The vibrational dynamics of solids is described by phonons constituting basic collective excitations in equilibrium crystals. Here, we consider a non-equilibrium active solid, formed by self-propelled particles, which bring the system into a non-equilibrium steady-state. We identify novel vibrational collective excitations of non-equilibrium (active) origin, which coexist with phonons and dominate over them when the system is far from equilibrium. These vibrational excitations are interpreted in the framework of non-equilibrium physics, in particular, stochastic thermodynamics. We call them “entropons” because they are the modes of spectral entropy production (at a given frequency and wave vector). The existence of entropons could be verified in future experiments on dense self-propelled colloidal Janus particles and granular active matter, as well as in living systems, such as dense cell monolayers.

## I. INTRODUCTION

Active matter^{1,2} includes a broad range of systems composed of particles that locally convert energy from the environment into directed motion.^{3,4} The energy exchange with the environment, often induced by chemical reactions or self-imposed gradients, leads to self-propulsion of the particles and internally drives an active system out of equilibrium.

Dense systems of self-propelled particles are rather ubiquitous in nature and often form crystalline structures. Cell monolayers of the human body^{5–8} and dense colonies of bacteria^{9–11} are common examples. Moreover, active colloidal Janus particles may cluster and form dense crystallites,^{12} named “living crystals.”^{13,14} These active crystals show fascinating phenomena uncommon for equilibrium solids^{15–18} ranging from intrinsic spatial velocity correlations,^{19–22} spontaneous velocity alignment,^{23–25} and non-Gaussian velocity distributions^{25,26} to “traveling” crystals,^{27–29} collective rotations,^{16,30–32} and even flocking.^{33} Activity also shifts the equilibrium freezing transition,^{34–37} affects the nature of the two-dimensional melting scenario,^{35,38–41} and changes the relaxation properties of amorphous solids.^{42–47} Systematic studies of collective excitations in active solids are still in their infancy,^{48} and they remain poorly explored even in general non-equilibrium solids (beyond active matter). Unveiling how the nature of the vibrations of a solid changes out of equilibrium represents an open issue relevant in statistical as well as solid-state physics.

External forces or internal mechanisms that dissipate energy drive a system away from equilibrium and spontaneously produce entropy. While self-propulsion is generated by an uptake of energy from the environment, likewise active particles dissipate energy, leading to local entropy production.^{49–54} Quantifying the non-equilibrium character of active systems via this observable has represented a topic of central interest in recent years.^{55–58} This subject has been investigated numerically, simulating both active field theories,^{59–62} active particle dynamics in interacting systems,^{63–65} as well as colloids in the presence of an active bath.^{66,67} Particular attention has been devoted to phase-separated configurations where the main contribution to the spatial profile of the entropy production has been observed at the interface between dense and dilute phases.^{59,68,69} Conversely, analytical results for entropy production have been only derived for simple cases, such as the potential-free particle,^{70–72} and for near-equilibrium regimes through perturbative methods.^{73,74} Entropy production in active crystals, however, remains unexplored.

In this paper, we fill this gap by discovering novel vibrational excitations characterizing non-equilibrium active solids and responsible for entropy production. Following the standard nomenclature of solid-state physics,^{75} we term these new modes “entropons” because each of them represents a spectral contribution to the entropy production rate. Unlike phonons, which describe the vibrational dynamics of equilibrium solids, entropons exist only in non-equilibrium, i.e., they are purely induced by activity. Entropons coexist with phonons but dominate over them for large activity and, therefore, represent the thermodynamically most relevant modes of an active crystal far from equilibrium. The underlying basic picture for active solids is shown in Fig. 1: in equilibrium, the thermal bath coupled to a crystal induces collective vibrational excitations like phonons (yellow color in Fig. 1); in non-equilibrium active solids, the self-propelling force injects energy into the crystal, this energy is dissipated in the environment, and the system produces entropy. In this process, entropons are generated as new collective vibrations, encoding the entropy production (orange color in Fig. 1). Our analysis is based on analytical theory combined with particle-resolved computer simulations and can, in principle, be verified in experiments with dense self-propelled Janus colloids^{12,76} or vibrated granular particles,^{29,77–79} as well as in living systems, such as confluent cell monolayers,^{6,8} which can exhibit crystalline patterns.

## II. MODEL

*N*inertial active Brownian particles (ABP) in a square box of size

*L*with periodic boundary conditions. Each particle with mass

*m*evolves through underdamped dynamics

^{80–83}for its position,

**x**

_{i}, and velocity $vi=x\u0307i$,

*ζ*_{i}and

*η*

_{i}are Gaussian white noises with zero average and unit variance. The term $fia=\gamma v0ni$ models the active force, with

*v*

_{0}being the swim velocity and

**n**

_{i}= (cos

*θ*

_{i}, sin

*θ*

_{i}) being the orientational unit vector, determined by an orientational angle

*θ*

_{i}. The coefficients

*γ*and

*T*correspond to the friction coefficient and temperature of the solvent bath, respectively, and define the inertial time

*τ*

_{I}=

*m*/

*γ*.

*D*

_{r}is the rotational diffusion coefficient, which determines the persistence time,

*τ*= 1/

*D*

_{r}, i.e., the time that the particle needs to randomize its orientation. The single-particle dynamics is often described in terms of the so-called active temperature, $Ta=v02\tau \gamma /2$, that vanishes in the equilibrium limits, either

*τ*→ 0 or

*v*

_{0}→ 0. The interaction force

**F**

_{i}stems from a soft repulsive pair potential,

*U*

_{tot}= ∑

_{i}

*U*(|

**x**

_{i}−

**x**

_{i}|), where $U=4\u03f5[(d0/r)12\u2212(d0/r)6]+\u03f5$ if

*r*<

*d*

_{0}2

^{1/6}and zero otherwise (Weeks–Chandler–Andersen potential

^{84}), with

*ϵ*and

*d*

_{0}being the energy scale and the particle diameter, respectively. The packing fraction $\varphi =\rho d02\pi /4=1.1$ (where

*ρ*=

*N*/

*L*

^{2}denotes the number density, with

*N*= 10

^{4}and

*L*/

*σ*= 88.6) is chosen high enough to ensure a solid-like behavior characterized by a defect-free triangular lattice, as illustrated in Fig. 1. Further details are reported in Appendix C. In the following, we chose an inertial regime such that

*τ*

_{I}∼

*τ*. Indeed, the inertial regime will play a key role in the spectral analysis of collective excitations.

## III. RESULTS

### A. Collective excitations and spectral entropy production

^{55,56,85}respectively (see Appendix A for details and Ref. 86 for a general review). The steady-state entropy production rate can be decomposed into its space-time Fourier spectrum as

**q**is the wave vector,

*ω*is the frequency, and Ω represents the area of the first two-dimensional Brillouin zone.

*σ*(

*ω*,

**q**) is analytically predicted as (see Appendix A)

*σ*as a function of the model parameters since the shape function

*K*(

*ω*) reads

*T*) and consists of a single term corresponding with phonons. In an active solid, $Cu\u0302u\u0302$ is formed by the sum of two contributions with thermal and active origins from which we can identify two distinct collective excitations. We have

*σ*(

*ω*,

**q**)/

*ω*

^{2}coincides with the spectral entropy production given by Eq. (3) and $Im[Ru\u0302u\u0302]=Im[Ru\u0302u\u0302(\omega ,q)]$ is the imaginary part of the response function due to a small perturbation,

*h*, evaluated in the frequency and wave vector domains. The response is defined as $Ru\u0302u\u0302(\omega ,q)=\delta \u27e8u\u0302(\omega )\u27e9/\delta h(\omega )|h=0$, and one has

*τ*(or

*v*

_{0}) at which the contribution of entropons becomes negligible compared to that of phonons. In the limit of vanishing active force (

*T*

_{a}→ 0), the response balances the lhs of Eq. (5) and the entropy production vanishes at equilibrium, as required. As a consequence, entropons disappear and thermal phonons remain the only collective excitations. Formula (5) can be interpreted as a Harada–Sasa relation applied to active solids at a given wave vector

^{87}(see also Ref. 59 for a derivation of the Harada–Sasa relation in active field theories) and justifies the decomposition into conventional phonons with thermal origin and entropons, arising from entropy production.

The decomposition of the collective excitations in phonons and entropons has a deep physical meaning directly linked to emergent collective phenomena. Indeed, active systems at high density show velocity patterns and spatial velocity correlations, independently observed in numerical simulations in Refs. 20 and 23 and in experiments based on cell monolayers in Ref. 8. Since thermally excited phonons do not produce spatial patterns in real space, entropons are clear evidence of novel excitations in active solids.

### B. Properties of entropons

A typical shape of *σ*(*ω*, **q**)*T*/*T*_{a} is shown in Fig. 2(a) as a function of *ω*/*ω*_{E} for a given **q**. A sharp peak occurs at a characteristic frequency *ω**(**q**). We identify this peak with an elementary excitation in the crystal and coin the term “entropon” to describe it, following the standard notation of elementary excitations in solids:^{75} Each entropon is identified with a peak of *σ*(*ω*, **q**).

Figures 2(b) and 2(c) show *σ*(*ω*, **q**)*T*/*T*_{a} as a function of *ω*/*ω*_{E} for different values of **q**, revealing a good agreement between theory, Eq. (3), and numerical simulations. Close to the equilibrium, in the regime of small persistence time such that *τ* = 1/*D*_{r} ≪ 1/*ω*_{E} [Fig. 2(b)], the peaks of *σ*(*ω*, **q**)*T*/*T*_{a} occur at the phonon frequency $\omega *(q)=\omega \u0304(q)$. In this regime, entropons have the same properties of phonons, but their amplitudes are small and proportional to *τ* (because of the prefactor *T*_{a}): The active force behaves as an additional thermal source at effective temperature *T*_{a}. In the opposite regime of large persistence time, *τ* = 1/*D*_{r} ≫ 1/*ω*_{E} [Fig. 2(c)], the peaks of *σ*(*ω*, **q**) are shifted to $\omega *(q)<\omega \u0304(q)$. As a consequence, entropons are different from phonons since the crystal vibrations are now peaked at frequencies not coinciding with those typical of equilibrium solids. The frequency *ω**(**q**) that maximizes *σ*(*ω*, **q**) is reported in Fig. 2(d) as a function of **q** for different rescaled persistence time, *τω*_{E}, while the difference $\omega \u0304(q)\u2212\omega *(q)$ is shown in Fig. 2(e) as a function of *τω*_{E} for several values of **q**. Despite *ω**(**q**) linearly increasing with **q** in the small persistence regime, a clear discrepancy from the linear law emerges in the large persistence regime for small **q**, where entropons follow a non-linear dispersion relation. As a consequence, the difference $\omega \u0304(q)\u2212\omega *(q)$ grows when *τω*_{E} is increased much more as **q** is decreased. Finally, the width of the peaks of *σ*(*ω*, **q**) increases with *τω*_{E} [compare the curves with the same colors in Figs. 2(b) and 2(c)], implying shorter-lived excitations at high activities.

*ω*of the spectral entropy production, $s(q)=\u222bd\omega \sigma (\omega ,q)=\rho TaT(\tau I+\tau )\u22121(1+\tau 2\tau I\tau +\tau I\omega \u03042(q))\u22121$, as a function of

*τω*

_{E}for different

**q**to quantify the weight of each entropon. For

*τω*

_{E}≲ 1, this observable is nearly

**q**-independent and increases with

*τω*

_{E}: The larger the

*τω*

_{E}, the larger is the contribution of each entropon to

*s*(

**q**). This linear behavior is due to the increase in the prefactor

*T*

_{a}∼

*τ*. For

*τω*

_{E}≳ 1, the value of $s(q)T/(v02\gamma )$ strongly depends on

**q**, displaying higher values for smaller

**q**. Interestingly, $s(q)T/(v02\gamma )$ decreases as

*τω*

_{E}is increased and, consequently, shows a non-monotonic behavior with a maximum that shifts for larger

*τω*

_{E}as

**q**is decreased. To shed light on this surprising non-monotonicity, the height of the peak of the rescaled spectral entropy production,

*σ*(

*ω**,

**q**)

*T*/

*T*

_{a}, is shown in Fig. 2(g) vs

*τω*

_{E}for different values of

**q**, revealing approximately the profile $\u223c1/(1+b\tau 2\omega E2)$, where

*b*=

*b*(

**q**) is a fitting parameter. Considering formula (3), the height of the peak of

*σ*is roughly determined by

*σ*(

*ω**,

**q**)

*T*/

*T*

_{a}∼

*K*(

*ω**). By approximating $\omega \u223c\omega \u0304(q)$, we obtain

*ξ*is the correlation length of the spatial velocity correlation, ⟨

**v**(

**r**) ·

**v**(0)⟩, of an active solid and reads

^{20}

**q**-dependence is negligible, we determine the typical wave vectors at which

*σT*/

*T*

_{a}starts decreasing,

**q**≲

**q**

_{*}∼ 1/

*ξ*(length scales larger than

*ξ*) provide the main contribution to the rescaled entropy, while those with

**q**≳

**q**

_{*}∼ 1/

*ξ*(length scales smaller than

*ξ*) have a much smaller weight. Since

*ξ*increases with

*τω*

_{E}, the modes with larger

**q**give a smaller contributions as

*τω*

_{E}is increased.

Finally, we recall that the present analysis has been obtained in the underdamped regime, when the inertial time *τ*_{I} is of the same order or larger than the typical relaxation time of the solid, given by the inverse of the Einstein frequency 1/*ω*_{E}. Indeed, in the overdamped regime for *τ*_{I} ≪ 1/*ω*_{E}, the spectrum shows only a single peak around *ω* ≈ 0, similarly to the case of equilibrium solids. This can be understood by looking at the denominator of Eq. (5) or Eq. (6) that suppresses the dependence from $\omega \u03042(q)$ for vanishing *τ*_{I}.

### C. Total entropy production rate

*ω*and

**q**and our prediction for

*σ*(

*ω*,

**q**) [Eq. (2)], we can derive analytically the expression for the global density of entropy production rate of the solid, $s\u0307$,

*T*) and is a function of persistence time

*τ*and inertial time

*τ*

_{I}. The integral in Eq. (10) can be analytically expressed in terms of elliptic functions depending on the parameters of the model (see Appendix B). Figure 3(a) plots $s\u0307$ as a function of

*τω*

_{E}, showing a non-monotic behavior: in the small persistence regime (

*τω*

_{E}≪ 1), $s\u0307\u223c0$ because the system behaves as an inertial solid in equilibrium. Increasing

*τω*

_{E}, the system departs from equilibrium and $s\u0307$ grows until reaching a maximum, roughly at

*τω*

_{E}∼ 1. For further increase in

*τω*

_{E}, $s\u0307$ decreases almost to zero. This decrease is consistent with the arrested states observed in dense active systems when

*τ*→

*∞*

^{44,88}for which $s\u0307\u22480$.

*ξ*. Here, we report the scaling behavior of $\sigma \u0307$ for

*ξ*≫ 1,

*ξ*increases suggests that the onset of spatial velocity correlations, characterizing active solids, reduces the entropy production rate. Coherent domains of strongly correlated velocities produce less entropy, while the incoherent behavior of the velocities of the particles determines a higher dissipation as if coherent motion, somehow, minimizes the effective friction between different particles. This disorder-order picture has been employed to explain non-monotonic behaviors also in conservative dynamics displaying collective motion.

^{89}

## IV. DISCUSSION

In conclusion, we have predicted new elementary vibrational excitations in active solids, termed entropons, because they are modes of the spectral entropy production. Our combined numerical and theoretical study revealed the properties of entropons, which coexist with phonons but dominate over phonons far from equilibrium when the active temperature is larger than the thermal one.

The concept of “entropons” as additional lattice vibrations has a broad generality that goes beyond monodisperse active crystals. It will certainly apply to binary crystals composed of active and passive particles.^{90,91} Moreover, we expect that in disordered dense systems, such as active glasses^{43,45–47,92,93} and dense active liquid crystals,^{94} entropons could play the dominant role of system excitations in determining entropy production. For instance, they could shed light on the activity-induced shift of the glass transition temperature.

Many experimental realizations of active crystals are available. Examples include confluent cell monolayers,^{6,8,95} dense assemblies of active colloids,^{76} as well as highly packed active granular systems,^{29,77,79} for which the solid structure has been achieved by connecting Hexbug particles by springs.^{48} Therefore, the existence of entropons can, in principle, be verified by analyzing particle trajectories in real space. To identify the contribution of entropons, it is crucial to have experimental access to the measurements of velocities and active forces, i.e., orientational angles. In this way, one can directly calculate the spectral entropy production in experiments by using Eq. (3). Alternatively, entropons can be measured indirectly without knowing the active forces through Eq. (5) by evaluating the dynamical correlations of the particle displacement and subtracting the contribution of thermal phonons, i.e., the response function due to a small perturbation.

Despite derived in the active case, entropons will characterize more general non-equilibrium crystals, as shown in Appendix E for a class of non-active solids driven out of equilibrium by temperature gradients. Therefore, the present paper reveals how concepts in the field of stochastic thermodynamics could have a determinant role in solid physics. For instance, we strongly believe that our theoretical approach can be applied to active solids characterized by polar, nematic, or even implicit alignment interactions between velocity and self-propelled angle. As a consequence, we speculate that entropons could be helpful to interpret and systematically analyze the collective excitations in a solid consisting of active granular particles, as the one experimentally realized in Ref. 48. Future studies on entropons can be manifold, for instance, the scattering of entropons near crystalline defects and the shifted spectra of entropons in a crystal exposed to a temperature gradient. This could lead to possible applications such as shock absorbers and active mass dampers as well as controlled heat radiators obtained by active crystals.

## ACKNOWLEDGMENTS

L.C. acknowledges support from the Alexander Von Humboldt foundation. U.M.B.M. and A.P. acknowledge support from MIUR PRIN 2017 Project No. 201798CZLJ. H.L. acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through the SPP 2265 under Grant No. LO 418/25-1.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Lorenzo Caprini**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). **U. Marini Bettolo Marconi**: Formal analysis (supporting); Investigation (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). **A. Puglisi**: Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). **Hartmut Löwen**: Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: DERIVATION OF THE SPECTRAL ENTROPY PRODUCTION

#### 1. Dynamics in the Fourier space

**v**

_{i}, and the active force $fia$,

*t*is the time-window numerically used to define the Fourier transform, corresponding to the time window of the simulations.

#### 2. Definition and implicit expression for the spectral entropy production

*σ*(

*ω*,

**q**) can be operatively calculated by using path-integral methods in frequency and wave vector domains starting from the path-integral definition of the entropy production rate $s\u0307$,

^{49,55,56,86}

*P*and

*P*

_{r}are the forward and backward trajectory expressed in terms of the forward and backward actions,

*A*and

*A*

_{r}, respectively, as

*L*(

*ω*,

**q**) and

*L*

_{r}(

*ω*,

**q**) are related to the forward and backward systems in

**q**and

*ω*, respectively. The former reads

**u**(

*ω*,

**q**) only. In this way,

*L*

_{r}(

*ω*,

**q**) is given by

*t*, the total entropy production rate $s\u0307$ can be expressed as follows (after some algebraic manipulations):

*t*and taking the limit

*t*→

*∞*, both terms vanish. However, the dependence on the force is implicitly contained in the expression for the entropy production rate, for instance, in the expression for the Einstein frequency [see Eq. (A15)].

##### a. Discretization scheme for the path-integral approach

To employ the path-integral approach, i.e., in writing the expressions for $A$ and $Ar$, we have adopted a continuous time formalism, assuming the Itô convention. Of course, since our system does not involve any multiplicative noise, the resulting entropy production rate is not affected by the choice of the integration scheme: Itô and Stratonovich conventions lead to the same results. The difference between Itô and Stratonovich conventions in the calculation of entropy production rate using path-integral approaches is described in detail in Refs. 49 and 85.

#### 3. Explicit expression for the spectral entropy production

To obtain an analytical expression for *σ*(*ω*, **q**), it is necessary to calculate the cross-correlation involved in Eq. (A11). In order to do that, we shall make two simplifying assumptions in the dynamics, Eq. (1) of the main text:

- Each particle performs small oscillations around a node of a triangular lattice so that the total inter-particle potential can be approximated as the sum of quadratic terms. Introducing the displacement
**u**_{i}of the particle*i*with respect to its equilibrium position, $xi0$, namely, $ui=xi\u2212xi0$, the pair potential, in the harmonic approximation, readswhere(A14)$Utot\u2248m\omega E22\u2211i\u2260j(uj\u2212ui)2,$*ω*_{E}is the Einstein frequency of the solid, given by(A15)$\omega E2=12mU\u2033(x\u0304)+U\u2032(x\u0304)x\u0304.$*ω*_{E}is determined by the first and second derivative of the interacting potential,*U*, calculated at $x\u0304$, i.e., the lattice constant of the solid, e.g., the average distance between neighboring particles in the solid. - We approximate the active Brownian particle (ABP) active force, $fia$, as an Ornstein–Uhlenbeck process for each particle evolving aswhere(A16)$\tau f\u0307ia=\u2212fia+\gamma v02\tau \zeta i,$
*ζ*_{i}is a vector of white noises such that ⟨*ζ*_{i}(*t*)*ζ*_{i}(0)⟩ =*δ*(*t*).

*d*/

*dt*→ −

*iω*, it is possible to rewrite the dynamics in a compact form upon introducing the vector $a\u0302(\omega ,q)=(u\u0302(\omega ,q),v\u0302(\omega ,q),f\u0302a(\omega ,q))$,

*i*being the imaginary unit and $I$ being the identity matrix.

**M**is the dynamical matrix that rules the deterministic part of the dynamics that reads

*T*stems for transpose matrix and the noise matrix

**reads**

*b*^{96–98}

*σ*(

*ω*,

**q**), i.e., Eq. (A11),

*τ*

_{I}=

*m*/

*γ*, the active temperature $Ta=\gamma v02\tau $, and the shake function

*K*(

*ω*) = 1/(1 +

*τ*

^{2}

*ω*

^{2}). From this equation, we can straightforwardly identify the spectral entropy production as

#### 4. Derivation of the positional dynamical correlation

*T*) and is activity-independent, while the second term has a pure active origin being proportional to the active temperature.

*σ*(

*ω*,

**q**)/

*ω*

^{2}] after considering the shape function

*K*(

*ω*) as in Eq. (4) of the main text. The first term in the right-hand side of Eq. (A30) is proportional to the imaginary part of the response function, $Ru\u0302u\u0302$, obtained by perturbing the system by a small force

*h*and defined as

^{87}i.e., a formula constraining dissipation (through entropy production), to response function and the correlation function of the displacement. At variance with the Harada–Sasa relation, developed in the case of particles in contact with a heat bath and driven out of equilibrium by ratcheting or constant forces, we have applied such a relation to active solids both in frequency and wave vector domains, and we have used such a relation to predict the existence of additional vibrational excitations, coexisting with phonons.

### APPENDIX B: DERIVATION OF THE TOTAL ENTROPY PRODUCTION OF THE SYSTEM

The prediction for the spectral entropy production *σ*(*ω*, **q**) allows us to calculate the entropy production of the system, $s\u0307$, explicitly as a function of the parameters of the model.

*r*(

**k**) for a triangular lattice in two-dimensions reads

**k**are given in terms of the cartesian components of

**q**by the relations

*z*in Eq. (B2) is given by

*K*(

*a*) is the modified Bessel function of the first kind. The term $s\u0307free$ coincides with the entropy production rate of a potential-free particle and is given by

*ξ*is the correlation length of the spatial velocity correlation given by

*z*,

*c*, and

*a*are explicit functions of the parameters of the model through the expression for

*ξ*and are given by

*ξ*, one gets Eq. (11) of the main text.

### APPENDIX C: NUMERICAL DETAILS

#### 1. Numerical details of the simulations

*δt*= 10

^{−5}

*τ*, that explicitly reads

**W**

_{i}and

**Y**

_{i}are Gaussian random numbers with unit variance and zero average. Simulations were usually performed up to a final time of at least

*t*

_{max}= 10

^{3}

*τ*. The active solid contains

*N*= 10

^{4}particles, at packing fraction

*ϕ*= 1.1, so that the size of the box reads

*L*= 85

*d*

_{0}. Positions are rescaled by the nominal particle diameter

*d*

_{0}, while time is rescaled by using Einstein frequency,

*ω*

_{E}(see definition in the main text). The dynamics is characterized by three typical times: the persistence time

*τ*(

*τ*= 1/

*D*

_{r}), the inertial time

*τ*

_{I}=

*m*/

*γ*, and the inverse of the Einstein frequency 1/

*ω*

_{E}(typical of equilibrium solids). Their combinations give rise to two dimensionless parameters:

*τω*

_{E}and

*τ*

_{I}

*ω*

_{E}. In the main text, the latter is kept fixed to

*τ*

_{I}

*ω*

_{E}= 3.5, and we have evaluated only the effect of the former parameter because it controls the non-equilibrium effects related to the spectral entropy production. In addition, the dynamics is characterized by the ratio between thermal and active temperature,

*T*/

*T*

_{a}= 10

^{−4}, where $Ta=v02\gamma \tau $ (as introduced and discussed in the text). Finally, the rescaled swim velocity reads

*v*

_{0}/(

*d*

_{0}

*ω*

_{E}) = 2 × 10

^{−1}, and the rescaled energy constant is given by $\u03f5/(md02\omega E2)=10\u22122$.

#### 2. Measure of the spectral entropy production and error estimate

To measure the spectral entropy production *σ*(*ω*, **q**), we have considered Eq. (3) (first line) of the main text. To do so, we have calculated the Fourier transform of data both in frequency *ω* and wave vector **q**, applying the definition on displacement, velocity, and active force by using Eqs. (A1a)–(A1c). Then, we have calculated the average dynamical correlation in the steady-state performing the average over the realization of the noise.

In the main text, Figs. 5(b) and 5(c) plot the spectral entropy production *σ*(*ω*, **q**) as a function of *ω*/*ω*_{E} for three values of the rescaled wave vector *qd*_{0} and two values of the rescaled persistence time *τω*_{E}. Data are reported without error bars for presentation reasons. Figure 2 contains the same data reported in Fig. 5 of the main text but split into several panels to confirm the agreement within the statistical error between data (points with error bars) and theoretical predictions (solid lines), given by Eq. (3) of the main text. In particular, panels (a)–(c) are realized with *τω*_{E} = 7 [Fig. 5(c) of the main text], while panels (d)–(f) are realized with *τω*_{E} = 7 × 10^{−2} [Fig. 5(b) of the main text]. Finally, panels (a) and (d), (b) and (e), and (c) and (f) plot *qd*_{0} = 0.1, 0.2, and 0.5, respectively.

### APPENDIX D: VISUALIZATION OF THE SYSTEM: SNAPSHOTS AND MEAN-SQUARE DISPLACEMENT

In this section, we provide a visualization of the system in real space. In Figs. 3(a)–3(d), we report several snapshot configurations in the plane of motion of the system for different values of the reduced persistent time *τω*_{E} = 7 × 10^{1}, 7, 7 × 10^{−1}, 7 × 10^{−2}. All of them display the typical hexagonal order characterizing two-dimensional crystals, as explicitly illustrated in Fig. 3, while the color gradient confirms the absence of spatial order in the direction of the active force, as expected.

*t*), defined as

*t*) is plotted as a function of

*tω*

_{E}for the same values of

*τω*

_{E}shown in the snapshot configurations (same column). After transient regimes, for instance, showing the typical subdiffusive effects characterizing the dynamics of solids, the system reaches a diffusive regime ∼

*t*outlined by the solid black line (which is a guide for the eyes). This analysis also confirms that the time window considered for the statistical analysis of the entropy production is enough to reach the steady state.

In Fig. 6, we also report the displacement maps Δ**x**_{i}(*t*) = **x**_{i}(*t*) − **x**_{i}(0) for different values of the reduced persistent time *τω*_{E} and for different duration times Δ*tω*_{E}, where Δ*t* = *t* − *t*_{0}. For the smaller value of *τω*_{E}, the displacement maps show weak structures that become more evident when Δ*tω*_{E} increases. Spatial structures in the displacement maps are more evident and larger when *τω*_{E} increases.

### APPENDIX E: THE CONCEPT OF ENTROPONS BEYOND ACTIVE SOLIDS

In this section, we show the generality of the concept of entropons that, in general, characterize a broad class of far from equilibrium solids, governed by a breaking of the fluctuation–dissipation theorem. This emphasizes the generality of entropons beyond active systems supporting our claim of generality expressed in the last sentence of the main paper.

*x*and

*y*of the displacement are driven by thermal baths at temperatures

*T*

_{x}and

*T*

_{y}, respectively (with

*T*

_{x}≠

*T*

_{y}), and mutually interact through a linear force,

*x*

_{ij}= (

*x*

_{i}−

*x*

_{j}) and Δ

*y*

_{ij}=

*y*

_{i}−

*y*

_{j}. While the force

**F**

_{i}is symmetric for

*x*→

*y*and

*y*→

*x*, the matrix of the coupling constants is not symmetric for

*K*→

*λ*and

*λ*→

*K*.

*K*and

*λ*have a different physical meaning. Indeed,

*K*is the diagonal coupling, providing the amplitude of the usual vectorial harmonic force of components (

*x*

_{i}−

*x*

_{j},

*y*

_{i}−

*y*

_{j}). Physically, this term is a harmonic force that keeps fixed the solid structure. Instead, the constant

*λ*determines the amplitude of an additional force, usually absent in equilibrium solids, that has components (

*y*

_{i}−

*y*

_{j},

*x*

_{i}−

*x*

_{j}) and that couples the dynamics of the

*x*component to that of the

*y*component of the crystal. The force acting on the

*x*component is determined by the

*y*component of the displacement (and vice versa).

*K*represents the elastic constant of the solid, and

*λ*is a coupling constant mixing

*x*and

*y*components. The solid is pushed far from equilibrium by the temperature difference and the coupling between

*x*and

*y*coordinates, as known even in the case of two particles coupled by a harmonic spring.

^{99}

**q**, and frequency domain,

*ω*. Using definitions (A1a) and (A1b), we define the vector of displacement, $u\u0302(\omega ,q)=(X\u0302(\omega ,q),Y\u0302(\omega ,q))$, around the lattice positions and velocity $v\u0302(\omega ,q)=(v\u0302x(\omega ,q),v\u0302y(\omega ,q))$. The transformed dynamics (E2) in Fourier space read

*ω*(

**q**) is given by

*λ*= 0 but also if

*T*

_{x}=

*T*

_{y}(equilibrium limit).

*σ*(

*ω*,

**q**) is the spectral entropy production given by Eq. (E8), while $Im(RX\u0302X\u0302(\omega ,q))$ and $Im(RY\u0302Y\u0302(\omega ,q))$ are the imaginary parts of the two diagonal elements of the response matrix that explicitly read

The expression for the dynamical correlations, contained in Eq. (E11), provides the information on the vibrational excitations of the non-equilibrium solids described by Eq. (E2). This observable is formed by the sum of two terms: (i) an equilibrium term, proportional to the response matrix, that provides the contribution of phonons along *x* and *y* directions and (ii) a second term that is proportional to the entropy production of the system and can be identified as the contribution of entropons, being proportional to the spectral entropy production. Consistently with our picture, this term disappears in the equilibrium limit for *λ* = 0 and/or *T*_{x} = *T*_{y}, i.e., when the entropy production rate vanishes.

The results of this section show that the picture of entropons goes beyond the case of active solids and characterizes more, in general, non-equilibrium solids reaching a steady state. In addition, in the case of solids driven out of equilibrium by the presence of different thermal baths, new vibrational excitations, encoding the spectral entropy production, coexist with thermal phonons.

## REFERENCES

*Elementary Excitations in Solids Lectures on Protons, Electrons, and Plasmons*

*Principles of Condensed Matter Physics*

*Theory of Simple Liquids: With Applications to Soft Matter*

*Computer Simulation of Liquids*