Entropons as collective excitations in active solids

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, thus, drives an active system intrinsically out of equilibrium.

Dense systems of self-propelled particles are rather ubiquitous in nature and often form crystalline structures. Cells monolayers of the human body [5, 6, 7, 8] and dense colonies of bacteria [9, 10, 11] which populate the human skin 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, 16, 17, 18] ranging from “traveling” crystals [19, 20, 21], intrinsic velocity correlations [22, 23, 24, 25] and spontaneous velocity alignment [26, 27] to collective rotations [16, 28]. Activity also shifts the equilibrium freezing transition significantly [29, 30, 31, 32] and affects the nature of the two-dimensional melting transition [30, 33, 34, 35, 36].

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 back into the environment. The resulting energy conversion leads to local entropy production [37, 38, 39, 40, 41]. Quantifying the non-equilibrium character of active systems via entropy production has represented a topic of central interest in recent years [42, 43, 44, 45]. This topic has been investigated numerically, simulating both active field theories [46, 47, 48, 49], active particle dynamics in interacting systems [50, 51, 52], as well as colloids in the presence of an active bath [53, 54]. 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 [46, 55, 56]. Conversely, analytical results for entropy production have been only derived for simple cases, such as the potential-free particle [57, 58, 59], and for near-equilibrium regimes through perturbative methods [60, 61]. Entropy production in active crystals, however, remains unexplored.

In this Letter, we fill this gap and calculate the entropy production in an active crystal. Performing a systematic spectral analysis, we show the existence of novel collective vibrational excitations in non-equilibrium active solids which are responsible for entropy production. Following the standard nomenclature of solid-state physics [62], we term these new modes “entropons”. 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 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, the active force, acting on each particle and producing their self-propulsion, 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 highly packed self-propelled colloidal Janus-particles [63] or vibrated granular grains [64, 65, 66], as well as in living systems such as confluent cell monolayers [6, 8].

We study a two-dimensional crystal of $N$ inertial active Brownian particles, in a square box of size $L$ with periodic boundary conditions. Each particle with mass $m$ evolves through an underdamped dynamics [67, 68, 69, 70] for its position, $\mathbf{x}_{i}$ and velocity $\mathbf{v}_{i}=\dot{\mathbf{x}}_{i}$

where $\boldsymbol{\zeta}_{i}$ and $\eta_{i}$ are Gaussian white noises with zero average and unit variance. The term $\mathbf{f}^{a}_{i}=\gamma v_{0}\mathbf{n}_{i}$ models the active force, $v_{0}$ being the swim velocity and $\mathbf{n}_{i}=(\cos\theta_{i},\sin\theta_{i})$ the orientational unit vector, determined by an orientational angle $\theta_{i}$. The coefficients $\gamma$ and $T$ correspond to the friction coefficient and the temperature of the solvent bath, respectively, and define the inertial time $\tau_{I}=m/\gamma$. $D_{r}$ is the rotational diffusion coefficient which determines the persistence time, $\tau=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, $T_{a}=v_{0}^{2}\tau\gamma$ that vanishes in the equilibrium limits, either $\tau\to 0$ or $v_{0}\to 0$. The interaction force $\mathbf{F}_{i}$ stems from a soft repulsive pair potential, $U_{\text{tot}}=\sum_{i}U(|\mathbf{x}_{i}-\mathbf{x}_{i}|)$, where $U=4\epsilon[(d_{0}/r)^{12}-(d_{0}/r)^{6}]$ if $r<2^{1/6}$ and zero otherwise (WCA potential), $\epsilon$ and $d_{0}$ being the energy scale and the particle diameter, respectively. The packing fraction $\phi=\rho d_{0}^{2}\pi/4$ (where $\rho=N/L^{2}$ denotes the number density), 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 the Supplemental Material (SM).

The non-equilibrium properties of the system are investigated by applying path-integral techniques to calculate the total entropy production, $\dot{s}=\lim_{t\to\infty}\langle\log\left(\mathcal{P}/\mathcal{P}_{r}\right)\rangle/t$, where $\mathcal{P}$ and $\mathcal{P}_{r}$ represent the path probabilities of the forward and backward trajectory [42, 43], respectively, see SM for details and Ref. [71] for a general review. The steady-state entropy production can be decomposed into its space-time Fourier spectrum as

where $\mathbf{q}$ is the wave vector, $\omega$ the frequency and $\Omega$ represents the area of the first two-dimensional Brillouin zone.

As a main result of this Letter, $\sigma(\omega,\mathbf{q})$ is analytically predicted as (see SM)

where the symbol CC stands for complex conjugate. The vectors $\tilde{\mathbf{f}}_{a}(-\omega,-\mathbf{q})$ and $\tilde{\mathbf{v}}(\omega,\mathbf{q})$ are the Fourier transforms of active force and velocity fields in the frequency and wave vector domains (see SM for their definitions). The second line of Eq. (Entropons as collective excitations in active solids) is obtained in the limit of a harmonic crystal and expresses $\sigma$ as a function of the parameters of the model since the shape function $K(\omega)$ explicitly reads

The term $\omega^{2}(\mathbf{q})$ in Eq. (Entropons as collective excitations in active solids) denotes the phonon dispersion relation of equilibrium solids, whose expression is reported in the SM for a triangular lattice. In general, $\omega^{2}(\mathbf{q})\propto\omega_{E}$ where $\omega_{E}=\frac{1}{2m}\left(U^{\prime\prime}(\bar{x})+\frac{U^{\prime}(\bar{x})}{\bar{x}}\right)$ is the Einstein frequency of the solid determined by the derivative of the potential calculated at the average distance between neighboring particles, $\bar{x}\sim 1/\sqrt{\rho}$.

The dynamical correlation function, $\left\langle{\tilde{\mathbf{u}}}(\omega,\mathbf{q})\cdot{\tilde{\mathbf{u}}}(-\omega,-\mathbf{q})\right\rangle$, where $\tilde{\mathbf{u}}$ is the Fourier transform of the particle displacement with respect to the unperturbed position of its lattice, can be expressed in terms of $\sigma$ as (see SM)

where $\text{Im}[\mathcal{R}_{uu}(\omega,\mathbf{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 $\mathcal{R}_{uu}(\omega,\mathbf{q})=\delta\langle\hat{\mathbf{u}}(\omega)\rangle/\delta h(\omega)|_{h=0}$ and one has

The relation (5) indicates that $\left\langle{\tilde{\mathbf{u}}}(\omega,\mathbf{q})\cdot{\tilde{\mathbf{u}}}(-\omega,-\mathbf{q})\right\rangle$ is the sum of two terms corresponding to: i) thermally excited crystal vibrations, independent of activity that we identify with phonons and ii) additional vibrational contributions of the solid associated with the entropy production that are purely induced by the active force. We remark that the latter are dominant if $T_{a}=v_{0}^{2}\tau\gamma\gg T$, i.e. far from equilibrium: There exists a typical $\tau$ (or $v_{0}$) at which their contribution becomes negligible with respect to the one of the phonons. In the limit of zero active force ($T_{a}\to 0$), the response balances the l.h.s of Eq. (5) and the entropy production vanishes as equilibrium imposes.

A typical shape of $\sigma(\omega,\mathbf{q})T/T_{a}$ is shown in Fig. 2 (a) as a function of $\omega/\omega_{E}$ for a given $\mathbf{q}$. A sharp peak occurs at a characteristic frequency $\omega^{*}(\mathbf{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 [62]: Each entropon is identified with a peak of ${\sigma}(\omega,\mathbf{q})$.

Fig. 2 (b) and (c) show $\sigma(\omega,\mathbf{q})T/T_{a}$ as a function of $\omega/\omega_{E}$ for different values of $\mathbf{q}$ revealing a good agreement between numerical simulations and theory, Eq.(Entropons as collective excitations in active solids). Close to the equilibrium, in the regime of small persistence time such that $\tau=1/D_{r}\ll 1/\omega_{E}$ (Fig. 2 (b)), the peaks of $\sigma(\omega,\mathbf{q})T/T_{a}$ occur at the phonon frequency $\omega^{*}(\mathbf{q})=\omega(\mathbf{q})$. In this regime, entropons have the same properties of phonons but their amplitudes are small and proportional to $\tau$ (because of the prefactor $T_{a}$). In these conditions, the active force behaves as an additional thermal source at effective temperature $T_{a}$. In the opposite regime of large persistence time, $\tau=1/D_{r}\gg 1/\omega_{E}$, (Fig. 2 (c)), the peaks of $\sigma(\omega,\mathbf{q})$ are shifted to $\omega^{*}(\mathbf{q})<\omega(\mathbf{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 $\omega^{*}(\mathbf{q})$ which maximizes $\sigma(\omega,\mathbf{q})$ is reported in Fig. 2 (d) as a function of $\mathbf{q}$ for different values of the rescaled persistence time, $\tau\omega_{E}$, while the difference $\omega(\mathbf{q})-\omega^{*}(\mathbf{q})$ is shown in Fig. 2 (e) as a function of $\tau\omega_{E}$ for different values of $\mathbf{q}$. Despite $\omega^{*}(\mathbf{q})$ linearly increases with $\mathbf{q}$ in the small persistence regime, a clear discrepancy from the linear law emerges in the large persistence regime for small $\mathbf{q}$, where entropons follow a non-linear dispersion relation. As a consequence, the difference $\omega(\mathbf{q})-\omega^{*}(\mathbf{q})$ grows when $\tau\omega_{E}$ is increased much more as $\mathbf{q}$ is decreased.

Fig. 2 (f) displays the integral over $\omega$ of the spectral entropy production, $s(\mathbf{q})=\int d\omega\sigma(\omega,\mathbf{q})=\rho\frac{T_{a}}{T}(\tau_{I}+\tau)^{-1}(1+\frac{\tau^{2}\tau_{I}}{\tau+\tau_{I}}\omega^{2}(\mathbf{q}))^{-1}$, as a function of $\tau\omega_{E}$ for different $\mathbf{q}$ to quantify the weight of each entropons. This observable is nearly $\mathbf{q}$-independent for $\tau\omega_{E}\lesssim 1$ but increases with $\tau\omega_{E}$: The larger $\tau\omega_{E}$ the larger is the contribution of each entropon to the entropy production. This linear behavior is due to the increase of the prefactor $T_{a}\sim\tau$. For $\tau\omega_{E}\gtrsim 1$, the value of $s(\mathbf{q})T/(v_{0}^{2}\gamma)$ strongly depends on $\mathbf{q}$ displaying higher values for smaller $\mathbf{q}$. In addition, $s(\mathbf{q})T/(v_{0}^{2}\gamma)$ decreases as $\tau\omega_{E}$ is increased and, consequently, shows a non-monotonic behavior with a maximum that shifts for larger $\tau\omega_{E}$ as $\mathbf{q}$ is decreased.

To shed light on this non-monotonicity, the height of the peak of the rescaled spectral entropy production, $\sigma(\omega^{*},\mathbf{q})T/T_{a}$, is shown in Fig. 2 (g) versus $\tau\omega_{E}$ for different values of $\mathbf{q}$, revealing approximately the profile $\sim 1/(1+b\tau^{2}\omega_{E}^{2})$ where $b=b(\mathbf{q})$ is a fitting parameter. The main contribution to $\sigma(\omega^{*},\mathbf{q})T/T_{a}$ is due to frequencies lower than $\omega^{*}(\mathbf{q})$. Considering formula (Entropons as collective excitations in active solids), the heigth of the peak of $\sigma$ is roughly determined by $\sigma(\omega^{*},\mathbf{q})T/T_{a}\sim K(\omega^{*})$. By approximating $\omega\sim\omega(\mathbf{q})$, we obtain

where $\xi$ is the correlation length of the spatial velocity correlation, $\langle\mathbf{v}(\mathbf{r})\cdot\mathbf{v}(0)\rangle$, of an active solid and reads [23]

Evaluating the denominator of Eq. (7), and requiring that the $\mathbf{q}$-dependence is negligible, we determine the typical wave vectors at which $\sigma T/T_{a}$ starts decreasing

The wave-vectors with $\mathbf{q}\lesssim\mathbf{q}_{*}\sim 1/\xi$ provide the main contribution to the rescaled entropy, while those with $\mathbf{q}\gtrsim\mathbf{q}_{*}\sim 1/\xi$ have a much smaller weight. Since $\xi$ increases with $\tau\omega_{E}$, the modes with larger $\mathbf{q}$ give a smaller contributions as $\tau\omega_{E}$ is increased.

By integrating our prediction for $\sigma(\omega,\mathbf{q})$ (Eq. (2)) over $\omega$ and $\mathbf{q}$, we can derive analytically the expression for the global density of entropy production of the solid, $\dot{s}$.

where $\dot{s}_{\text{free}}=\rho\frac{T_{a}}{T}(\tau_{I}+\tau)^{-1}$ is the density of entropy production of non-interacting active particles. $\dot{s}_{\text{free}}$ is proportional to the ratio between the active temperature ($T_{a}=v_{0}^{2}\gamma\tau$) and the thermal temperature ($T$), and is a function of the dynamical properties of the system, such as persistence time $\tau$ and inertial time $\tau_{I}$. The integral in Eq. (10) can be analytically expressed in terms of elliptic functions in terms of the parameters of the model (see SM). Figure 3 (a) plots $\dot{s}$ as a function of $\tau\omega_{E}$ showing a non-monotic behavior: in the small persistence regime ($\tau\omega_{E}\ll 1$), $\dot{s}\sim 0$ because the system behaves as an inertial solid in equilibrium. Increasing $\tau\omega_{E}$, the system departs from equilibrium and entropy production grows until reaching a maximum, roughly at $\tau\omega_{E}\sim 1$. For further increase of $\tau\omega_{E}$, the value of $\dot{s}$ decreases until almost vanishes.

While the increase of $\dot{s}$ is expected when the system departs from equilibrium and is well-explained by the increase (up to saturation) of the non-interacting entropy production $\dot{s}_{\text{free}}\sim\tau/(\tau+\tau_{I})$, the physical picture behind the decrease in the large persistence regime can be understood by expressing $\dot{s}$ as a function of $\xi$. Here, we report the scaling behavior of $\dot{\sigma}$ for $\xi\gg 1$

which is shown in Fig. 3 (b). The decrease of $\dot{s}$ as $\xi$ increases suggests that the onset of spatial velocity correlations characterizing the solid at higher densities reduces the entropy production. 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.

In conclusion, we have predicted new elementary vibrational excitations in active solids, termed entropons, which emerge from the spectral entropy production and are generated uniquely by activity. Our combined numerical and theoretical study revealed the properties of entropons and demonstrated that they dominate over phonons far from equilibrium.

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 [72]. Moreover, we expect that in disordered dense systems, such as active glasses [73, 74, 75, 76] and dense active liquid crystals [77], 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 biological tissues [78], confluent cell monolayers [6, 8], dense assemblies of active colloids [63] as well as highly packed active granular systems [64, 66], for which the solid structure has been recently achieved by connecting Hexbug particles by springs [79]. Therefore the existence of entropons can in principle be verified by analyzing particle trajectories in real space.

Future studies on entropons can be manifold by including the scattering of entropons near crystalline defects, the shifted spectra of entropons in a crystal exposed to a temperature gradient. This could lead to possible applications like shock absorbers and active mass dampers as well as controlled heat radiators obtained by active crystals.