Thermodynamics of emergent structure in active matter

TL;DR

This paper investigates the nonequilibrium thermodynamics of active matter systems, revealing how entropy production varies across phases and identifying thermodynamic signatures of phase transitions and micro-structural features.

Contribution

It introduces a fully under-damped model for active matter, uncovering hidden entropy productions and analyzing their relation to phase behavior and structural defects.

Findings

Discontinuous entropy production at phase transitions

Identification of dissipation rates for different phases

Entropy production linked to micro-structural features like defects

Abstract

Active matter is rapidly becoming a key paradigm of out-of-equilibrium soft matter exhibiting complex collective phenomena, yet the thermodynamics of such systems remain poorly understood. In this letter we study the nonequilbrium thermodynamics of large scale active systems capable of mobility-induced phase separation and polar alignment, using a fully under-damped model which exhibits hidden entropy productions not previously reported in the literature. We quantify steady state entropy production at each point in the phase diagram, revealing characteristic dissipation rates associated with the distinct phases and configurational structure. This reveals sharp discontinuities in the entropy production at phase transitions and facilitates identification of the thermodynamics of micro-features, such as defects in the emergent structure. The interpretation of the time reversal symmetry in the dynamics of the particles is found to be crucial.

Figures

Tables

Table 1. Table 1: Expected entropy production rate for a system of interacting ABPs at steady state, described by over-damped and under-damped dynamics, and for odd and even interpretation of the particles’ heading under TRS.

$⟨\Delta {\dot{S}}^{\text{tot}}⟩$	Over-damped	Under-damped
TRS-odd heading	${\sum }_{a=1}^N(\beta {\left(I{\gamma }_{\text{R}}\right)}^1⟨{𝒯}_a^2(𝐫,\theta )⟩+{\left(I{\gamma }_{\text{R}}\right)}^1⟨{\partial }_{{\theta }_a}{𝒯}_a(𝐫,\theta )⟩$ $+{\sum }_{j=1}^2\beta v_0⟨{\partial }_{r_a^j}{U}_a\left(𝐫\right){\hat{e}}^j\left({\theta }_a\right)⟩)$	${\sum }_{a=1}^N({\gamma }_{\text{R}}(\beta I⟨{\omega }_a^2⟩-1)$ $+{\sum }_{j=1}^2\gamma (\beta m{v}_0^2/2-\beta m⟨{\left(v_a^j\right)}^2⟩+1))$
TRS-even heading	${\sum }_{a=1}^N(\beta {\left(I{\gamma }_{\text{R}}\right)}^1⟨{𝒯}_a^2(𝐫,\theta )⟩+{\left(I{\gamma }_{\text{R}}\right)}^1⟨{\partial }_{{\theta }_a}{𝒯}_a(𝐫,\theta )⟩$ $+m\gamma \beta v_0^2-{\sum }_{j=1}^2\beta v_0⟨{\partial }_{r_a^j}{U}_a\left(𝐫\right){\hat{e}}^j\left({\theta }_a\right)⟩)$	${\sum }_{a=1}^N({\gamma }_{\text{R}}(\beta I⟨{\omega }_a^2⟩-1)$ $+{\sum }_{j=1}^2\gamma (\beta m⟨{\left(v_a^j\right)}^2⟩-1))$

Equations

$$\begin{array}{rl}d{r}_a^j=& (m\gamma {)}^{\text{shortminus}1}{\mathcal{F}}_a^j(\mathbf{r},\theta )dt+\sqrt{2/\beta m\gamma }d{W}_{r_a^j},\end{array}$$

$$\begin{array}{rl}d{\theta }_a=& (I{\gamma }_{\text{R}}{)}^{\text{shortminus}1}{\mathcal{T}}_a(\mathbf{r},\theta )dt+\sqrt{2/\beta I{\gamma }_{\text{R}}}d{W}_{{\theta }_a},\end{array}$$

$$\begin{array}{rl}d{r}_a^j=& v_a^{j}dt,\end{array}$$

$$\begin{array}{rl}d{v}_a^j=& -\gamma v_a^{j}dt+m^{\text{shortminus}1}{\mathcal{F}}_a^j(\mathbf{r},\theta )dt+\sqrt{2\gamma /\beta m}d{W}_{v_a^j},\end{array}$$

$$\begin{array}{rl}d{\theta }_a=& {\omega }_{a}dt,\end{array}$$

$$\begin{array}{rl}d{\omega }_a=& -{\gamma }_{\text{R}}{\omega }_{a}dt+I^{\text{shortminus}1}{\mathcal{T}}_a(\mathbf{r},\theta )dt+\sqrt{2{\gamma }_{\text{R}}/\beta I}d{W}_{{\omega }_a},\end{array}$$

$$⟨\Delta {\dot{S}}^{\text{tot}}⟩=m\beta \gamma v_0^2(1-{\gamma }^2\mathcal{G})$$

$$⟨\Delta {\dot{S}}^{\text{tot}}⟩=m\beta {\gamma }^3{v}_0^2\mathcal{G}$$

$$\mathcal{G}=\underset{t\to \infty }{\lim }2{\int }_0^{t}d{t}_1{\int }_0^{t}d{t}_2{e}^{\text{shortminus}\gamma (2t\text{shortminus}{t}_1\text{shortminus}{t}_2)}⟨{\hat{e}}^1(t_1){\hat{e}}^1(t_2)⟩.$$

$$⟨\Delta {\dot{S}}_a^{\text{med}}|\Omega ⟩$$

$$+{\sum }_{j=1}^2\gamma \left(\beta m⟨(v_a^j-v_0{\hat{e}}^j({\theta }_a){)}^2⟩-1\right),$$

$$⟨\Delta {\dot{S}}_a^{\text{med}}|\Omega ⟩$$

$$d{x}_a=A_{x_a}(\Omega (t),t)dt+B_{x_a}(\Omega (t),t)d{W}_{x_a}.$$

$$d{x}_a=A_{x_a}^{\text{REV}}(\Omega (t),t)dt+A_{x_a}^{\text{IR}}(\Omega (t),t)dt+B_{x_a}(\Omega (t),t)d{W}_{x_a}.$$

$$d\Delta S^{\text{med}}=\sum _{a=1}^N\sum _{x_a\in \{r_a^1,r_a^2,v_a^x,v_a^2,{\theta }_a,{\omega }_a\}}\frac{A_{x_a}^{\text{IR}}(\Omega (t))}{D_{x_a}}\circ d{x}_a-\frac{A_{x_a}^{\text{IR}}(\Omega (t))A_{x_a}^{\text{REV}}(\Omega (t))}{D_{x_a}}dt$$

$$\begin{array}{rl}d{r}_a^j=& v_0{\hat{e}}^j({\theta }_a)dt-\frac{1}{m\gamma }{\partial }_{r_a^j}{U}_a(\mathbf{r})dt+\sqrt{\frac{2}{\beta m\gamma }}d{W}_{r_a^j},\end{array}$$

$$\begin{array}{rl}d{\theta }_a=& \frac{1}{I{\gamma }_{\text{R}}}{\mathcal{T}}_a(\mathbf{r},\theta )dt+\sqrt{\frac{2}{\beta I{\gamma }_{\text{R}}}}d{W}_{{\theta }_a}.\end{array}$$

$$d\Delta S^{\text{med}}=\sum_{a=1}^{N}\Bigg{(}\beta{\cal T}_{a}(\mathbf{r},\theta)\circ d\theta_{a}-\sum_{j=1}^{2}\beta\bigg{(}\partial_{{r}_{a}^{j}}U_{a}(\mathbf{r})\circ d{r}_{a}^{j}+v_{0}\hat{e}^{j}(\theta_{a})\partial_{{r}_{a}^{j}}U_{a}(\mathbf{r})dt\bigg{)}\Bigg{)}.$$

$$\frac{d\langle\Delta S^{\text{med}}\rangle}{dt}=\sum_{a=1}^{N}\Bigg{(}\frac{\beta\langle{\cal T}_{a}^{2}(\mathbf{r},\theta)\rangle}{I\gamma_{\text{R}}}+\frac{\langle\partial_{\theta_{a}}{\cal T}_{a}(\mathbf{r},\theta)\rangle}{I\gamma_{\text{R}}}+\sum_{j=1}^{2}\beta v_{0}\left\langle\partial_{{r}_{a}^{j}}U_{a}(\mathbf{r})\hat{e}^{j}(\theta_{a})\right\rangle\Bigg{)}.$$

$$d\Delta S^{\text{med}}=\sum_{a=1}^{N}\Bigg{(}\beta{\cal T}_{a}(\mathbf{r},\theta)\circ d\theta_{a}+\sum_{a=1}^{N}\beta\Big{(}v_{0}\hat{e}^{j}(\theta_{a})-\partial_{{r}_{a}^{j}}U_{a}(\mathbf{r})\Big{)}dt\Bigg{)}.$$

$$\frac{d\langle\Delta S^{\text{med}}\rangle}{dt}=\sum_{a=1}^{N}\Bigg{(}\frac{\beta\langle{\cal T}_{a}^{2}(\mathbf{r},\theta)\rangle}{I\gamma_{\text{R}}}+\frac{\langle\partial_{\theta_{a}}{\cal T}_{a}(\mathbf{r},\theta)_{a}\rangle}{I\gamma_{\text{R}}}+\beta m\gamma v_{0}^{2}-\sum_{a=1}^{N}\beta v_{0}\left\langle\partial_{{r}_{a}^{j}}U_{a}(\mathbf{r})\hat{e}^{j}(\theta_{a})\right\rangle\Bigg{)}.$$

$$\begin{array}{rl}d{r}_a^j=& v_a^{j}dt,\end{array}$$

$$\begin{array}{rl}d{v}_a^j=& -\gamma v_a^{j}dt+\gamma v_0{\hat{e}}^j({\theta }_a)dt-\frac{1}{m}{\partial }_{r_a^j}{U}_a(\mathbf{r})dt+\sqrt{\frac{2\gamma }{\beta m}}d{W}_{v_a^j},\end{array}$$

$$\begin{array}{rl}d{\theta }_a=& {\omega }_{a}dt,\end{array}$$

$$\begin{array}{rl}d{\omega }_a=& -{\gamma }_{\text{R}}{\omega }_{a}dt+\frac{1}{I}{\mathcal{T}}_a(\mathbf{r},\theta )dt+\sqrt{\frac{2{\gamma }_{\text{R}}}{\beta I}}d{W}_{{\omega }_a}.\end{array}$$

$$d\Delta S^{\mathrm{med}}$$

$$-\sum _{i=1}^2\beta m{v}_j^i\circ d{v}_j^i+\beta v_i^j({\partial }_{r_j^i}{U}_j(\mathbf{r})-v_0{\hat{e}}^i({\theta }_j))dt$$

$$\frac{d\Delta S^{\text{med}}}{dt}=\sum_{a=1}^{N}\Bigg{(}\gamma_{\text{R}}\Big{(}\beta I\left\langle\omega_{a}^{2}\right\rangle-1\Big{)}+\sum_{j=1}^{2}\gamma\bigg{(}\beta m\left\langle({v}^{j}_{a})^{2}\right\rangle-1\bigg{)}\Bigg{)},$$

$$\frac{d\Delta S^{\text{med}}}{dt}=\sum_{a=1}^{N}\Bigg{(}\gamma_{\text{R}}\Big{(}\beta I\left\langle\omega_{a}^{2}\right\rangle-1\Big{)}-\sum_{j=1}^{2}m\gamma v_{0}\beta\langle\hat{e}^{j}(\theta_{a})v_{a}^{j}\rangle\Bigg{)}$$

$$\sum _{j,i}⟨{\hat{e}}^i({\theta }_j)v_j^i⟩$$

$$d\Delta S^{\mathrm{med}}$$

$$+\sum _{i=1}^2\beta (-m{v}_j^i+m{v}_0{\hat{e}}^i({\theta }_j))\circ d{v}_j^i+\beta (-v_i^j+v_0{\hat{e}}^i({\theta }_j)){\partial }_{r_j^i}{U}_j(\mathbf{r})dt$$

$$d\Delta S^{\text{med}}=\sum_{a=1}^{N}\Bigg{(}\gamma_{\text{R}}\Big{(}\beta I\left\langle\omega_{a}^{2}\right\rangle-1\Big{)}+\sum_{j=1}^{2}\gamma\bigg{(}\beta m\left\langle(\mathbf{v}^{j}_{a}-v_{0}\hat{e}^{j}(\theta_{a}))^{2}\right\rangle-1\bigg{)}\Bigg{)},$$

$$d⟨\Delta S^{\mathrm{med}}⟩$$

$$\frac{d\langle\Delta S^{\text{med}}\rangle}{dt}=\sum_{a=1}^{N}\Bigg{(}\gamma_{\text{R}}\Big{(}\beta I\left\langle\omega_{a}^{2}\right\rangle-1\Big{)}+\sum_{j=1}^{2}\gamma\left(\frac{\beta mv_{0}^{2}}{2}-\beta m\left\langle({v}^{j}_{a})^{2}\right\rangle+1\right)\Bigg{)}$$

$$d⟨\Delta S^{\mathrm{med}}⟩$$

Full text

Active matter consists of particles that can consume stored free energy reserves in order to self-propel, and as such are characteristically out-of-equilibrium Ramaswamy (2010); Marchetti et al. (2013); Bechinger et al. (2016); Ramaswamy (2017). Examples encompass a wide range of systems, including self-catalytic colloidal suspensions Bialké et al. (2015), swimming bacteria Czirók et al. (2001); Sokolov et al. (2009), migrating cells Szabo et al. (2006) and animal groups Parrish et al. (2002); Buhl et al. (2006); Cavagna et al. (2010). Self-propulsion, in combination with interactions amongst the particles, can give rise to non-trivial collective dynamics not observed in matter at thermal equilibrium, such as gathering, swarming and swirling Vicsek and Zafeiris (2012).

Widely used models of active particles include Active Brownian Particles (ABPs) Schimansky-Geier et al. (1995) and Active Ornstein-Uhlenbeck Particles (AOUPs) Fodor et al. (2016). Collective motion and kinetic phase transitions can be observed in such models, with the introduction of volume exclusion, e.g., between two-dimensional discs Marchetti et al. (2016). Indeed, systems of both ABPs and AOUPs have been shown to exhibit motility-induced phase separation (MIPS), where the particles arrange themselves into regions of high and low density Fily and Marchetti (2012); Redner et al. (2013a, b); Cates and Tailleur (2015); Solon et al. (2015); Marchetti et al. (2016); Fodor et al. (2016); Cugliandolo et al. (2017); Digregorio et al. (2018); Siebert et al. (2018). In addition, a recent study has shown that ABPs with alignment interactions can exhibit polar collective motion (or ‘flocking’) Martín-Gómez et al. (2018).

Determining the phase diagrams for such behavior has been an active area of research Fily and Marchetti (2012); Redner et al. (2013a, b); Kapfer and Krauth (2015); Cugliandolo et al. (2017); Martín-Gómez et al. (2018); Digregorio et al. (2018); Siebert et al. (2018), however, there has been less focus on the thermodynamics, especially on the nonequilibrium character of the different kinetic phases, despite some progress in related field-theoretic models Nardini et al. (2017). For molecular approaches, Fodor et al. Fodor et al. (2016) investigated the entropy production in a system of AOUPs with no alignment interactions, arguing that in a harmonic trap the dynamics respect detailed balance such that the system is in an effective equilibrium. Later, Mandal et al. Mandal et al. (2017) demonstrated that when a different definition of entropy production is used the nonequilibrium character can be recovered. Recently, Shankar et al. Shankar and Marchetti (2018) investigated the ‘hidden’ components of the entropy production observed only in under-damped descriptions of the particles’ translational dynamics, reporting a key dependence on the time reversal symmetry (TRS) interpretation of the self-propulsion force, mirroring the distinct approaches of Fodor and Mandal. However, the study only considered free, non-interacting particles.

In contrast, we consider a large system of ABPs interacting via volume exclusion as well as alignment enabling investigation of the entropy production associated with the emergent collective motion. We derive expressions for the entropy production for both over and under-damped models, under odd and even interpretations of the parity of the particles’ heading. This reveals an additional hidden component due to coarse grained rotational dynamics even with under-damped translational dynamics. Simulation of the under-damped model allows us not only to construct the phase diagram, but also quantify the steady state thermodynamics at each point in the space. Further, we are able to examine the spatial distributions of entropy production associated with the distinct phases alongside micro-features, such as defects, in the emergent structures.

We consider a system of $N$ two-dimensional, disc-shaped ABPs of radius $R$, mass $m$, moment of inertia $I$, self-propulsion speed $v_{0}$ and translational and rotational mobility coefficients $\gamma$ and $\gamma_{\text{R}}$. The position and heading of each particle $a$ are denoted as $\mathbf{r}_{a}{=}\{r_{a}^{1},r_{a}^{2}\}$ and $\theta_{a}$ respectively (variables without subscripts or superscripts are to be understood as the total set of such variables in the system, e.g. $\mathbf{r}{=}\{\mathbf{r}_{1},\ldots,\mathbf{r}_{N}\}$). The self-propulsion force is modeled as ${\cal P}(\theta_{a}){=}\{{\cal P}^{1}(\theta_{a}),{\cal P}^{2}(\theta_{a})\}{=}m\gamma v_{0}\hat{e}(\theta_{a})$, where \hat{e}(\theta_{a}){=}$$\{\hat{e}^{1}(\theta_{a}),\hat{e}^{2}(\theta_{a})\}$${=}\{\cos(\theta_{a}),\sin(\theta_{a})\}. Excluded volume interactions are modeled using a truncated and shifted Lennard-Jones potential $U(\mathbf{r}){=}\sum_{a}U_{a}(\mathbf{r})$ with $U_{a}(\mathbf{r}){=}\sum_{b\neq a}\epsilon[(2R/\mathbf{r}_{ab})^{12}{-}(4R/\mathbf{r}_{ab})^{6}]+\epsilon$ if $|\mathbf{r}_{ab}|{\leq}R$ and $U_{a}(\mathbf{r}){=}0$ if $|\mathbf{r}_{ab}|{>}R$, where $\epsilon$ is the depth of the potential well and $\mathbf{r}_{ab}{=}\mathbf{r}_{a}{-}\mathbf{r}_{b}$. Finally, informed by the Kuramoto model Acebrón et al. (2005), alignment interactions are modeled as ${\cal T}_{a}(\mathbf{r},\theta){=}{-}K\sum_{b\neq a}g(\mathbf{r}_{ab})\sin(\theta_{a}{-}\theta_{b})$ where $K$ is the coupling strength and $g(\mathbf{r}_{ab}){=}1$ if $|\mathbf{r}_{ab}|{\leq}2R$ and zero otherwise.

A minimal description of the system is given by the following (over-damped) stochastic differential equations (SDEs):

[TABLE]

where ${\cal F}^{j}_{a}(\mathbf{r},\theta){=}{\cal P}^{j}(\theta_{a}){-}\partial_{r^{j}_{a}}U_{a}(\mathbf{r})$, $i$ is the spatial dimension, $\beta$ is the inverse temperature (with units $k_{B}{=}1$) and $W_{r^{j}_{a}}$ and $W_{\theta_{a}}$ are uncorrelated Wiener processes, such that $\langle dW_{r^{j}_{a}}dW_{\theta_{a}}\rangle{=}0$, $\langle dW_{\theta_{a}}dW_{\theta_{b}}\rangle{=}\delta_{ab}dt$ and $\langle dW_{r^{j}_{a}}dW_{r^{k}_{b}}\rangle{=}\delta_{jk}\delta_{ab}dt$. A finer grained description of the dynamics includes the translational and rotational velocities $\mathbf{v}_{a}{=}\{v_{a}^{1},v_{a}^{2}\}$ and $\omega_{a}$ through the under-damped SDEs:

[TABLE]

where $W_{v_{a}^{j}}$ and $W_{\omega_{a}}$ are also independent Wiener processes.

To understand the thermodynamics of these models we may turn to the framework of stochastic thermodynamics Seifert (2008). A central quantity of interest is the steady dissipation, or entropy production, which in such a formalism can be interpreted as a measure of dynamical irreversibility. Taking $k_{B}{=}1$ and defining $\Omega{=}\{\mathbf{r},\mathbf{v},\theta,\omega\}$ (or $\Omega{=}\{\mathbf{r},\theta\}$ for the over-damped system) as the total state of the system, the entropy production of an individual realization $\vec{\Omega}{=}\{\Omega(t)|t{\in}[t_{0},\tau]\}$, over the interval $[t_{0},\tau]$, is given by $\Delta S^{\text{tot}}{=}\ln(P[\vec{\Omega}]/P^{\dagger}[\vec{\Omega}^{\dagger}])$ Seifert (2005, 2008). Here ${P}$ and ${P}^{\dagger}$ are the probability measures for the forward and time reversed dynamics, $\vec{\Omega}^{\dagger}{=}\{\Omega^{\dagger}(t)|t{\in}[t_{0},\tau]\}$ and $\Omega^{\dagger}(t){=}\mathbf{\varepsilon}\Omega(\tau{+}t_{0}{-}t)$ where $\mathbf{\varepsilon}$ is a time reversal operator Spinney and Ford (2012a). Consequently, the entropy production is equal to the log ratio of the likelihood of a given trajectory against its time reverse. For stationary, autonomous and time symmetric dynamics, e.g. ABPs in a steady state, ${P}^{\dagger}{=}{P}$.

The total entropy production, comprising the change in entropy of the system and the environment, obeys an integral fluctuation theorem $\langle\exp[-\Delta S^{\text{tot}}]\rangle{=}1$. Thus the strict inequality $\langle\Delta S^{\text{tot}}\rangle{\geq}0$ holds by Jensen’s inequality, characterizing the second law. For SDEs, expressions for the total entropy production can be found exactly given knowledge of the probability density functions over the variables Spinney and Ford (2012), whilst expressions for the environmental entropy production can be determined in terms of the trajectories only. In the steady state, however, the expected change of system entropy vanishes and the mean medium entropy production is equal to the mean total entropy production allowing empirical calculation of expectations without the need for solving the associated Fokker-Planck equation.

Utilizing the formalism in Spinney and Ford (2012) we derive the expected medium entropy production for the total system of ABPs, assuming a steady state, for both over-damped and under-damped dynamics. The results depend crucially on the operator $\mathbf{\varepsilon}$, with uncertainty in the literature as to the time reversal symmetry of the particles’ orientation Shankar and Marchetti (2018). These entropy productions, under both odd and even interpretations of $\theta$ are reported in Table 1 with details in the Supplemental Material (SM). These expressions are quite general, however, in the absence of alignment, external, and exclusion interactions, such that ${\cal T}_{a}(\mathbf{r},\theta){=}0$ and $\partial_{r_{a}^{j}}U(\mathbf{r}){=}0$, we can recover and generalize the results for free ABPs in Shankar and Marchetti (2018). The over-damped results follow directly, however, results using under-damped translational dynamics depend on the treatment of the rotational degrees of freedom. Generally the individual free particle entropy production takes the form

[TABLE]

for the odd interpretation of $\theta$ and

[TABLE]

for the even interpretations of $\theta$, where

[TABLE]

Shankar et al. Shankar and Marchetti (2018), whilst considering under-damped translational motion, utilize over-damped rotational motion with $2\langle\hat{e}^{1}(t_{1})\hat{e}^{1}(t_{2})\rangle{=}e^{(\beta I\gamma_{\text{R}})^{-1}|t_{1}\shortminus t_{2}|}$ and thus $\mathcal{G}_{\text{over}}{=}(\gamma(\gamma{+}(\beta I\gamma_{\text{R}})^{\shortminus 1}))^{\shortminus 1}$. However, a fully under-damped description of the free-particle dynamics yields $2\langle\hat{e}^{1}(t_{1})\hat{e}^{1}(t_{2})\rangle{=}e^{(\beta I\gamma^{2}_{R})^{\shortminus 1}(1\shortminus\gamma_{\text{R}}|t_{1}\shortminus t_{2}|\shortminus\exp[\shortminus\gamma_{\text{R}}|t_{1}\shortminus t_{2}|])}$. Except for specific choice of parameters (e.g., all free parameters set to $1$), the integral has no closed form solution, but strictly satisfies $\mathcal{G}_{\text{under}}\!\geq\!\mathcal{G}_{\text{over}}$ indicating an additional hidden component in the entropy production (see details in the SM).

This illustrates the well known property that contributions are lost through coarse graining procedures Esposito (2012). Such absent terms have been referred to as ‘anomalous’ Celani et al. (2012) or ‘hidden’ and have been previously implicated in heat transfer where under-damped models are crucial in order to observe physically plausible entropy productions Spinney and Ford (2012). However, here the results are particularly nuanced as translational entropy productions in the free particle results are hidden due to coarse-graining in the rotational degrees of freedom and the discrepancies are non-trivial. For instance, setting all free parameters to $1$ yields $\mathcal{G}_{\text{under}}{=}(2e-4)\mathcal{G}_{\text{over}}\!\simeq\!1.44\;\mathcal{G}_{\text{over}}$, with commensurate over and under estimates in the entropy production for odd and even $\theta$ respectively. In light of this, despite apparent additional complication, we proceed utilizing the fully under-damped model so that we can be assured no features are either missing or are introduced as artifacts.

Importantly, as the Wiener processes are assumed to be uncorrelated (also known as a bipartite, or rather multipartite property  Spinney et al. (2018)), we may associate entropy productions with individual particles, with the total being their sum. For the odd interpretation the expected entropy production rate in the medium for particle $a$ (no longer assuming a steady state) is

[TABLE]

where again $k_{B}{=}1$. For the even interpretation we have

[TABLE]

The distinction between the parity interpretations is striking: under the even parity interpretation, the entropy production is manifestly a measure of the deviation away from equipartition expected at thermodynamic equilibrium in both the translational and rotational degrees of freedom. In contrast, under the odd parity interpretation the entropy production arising from the translational variables is modified such that it quantifies deviation from an effective equipartition, relative to the instantaneous heading and typical speed.

The system is simulated by integrating Eqs. (3-6) with a stochastic velocity Verlet algorithm (details can be found in the SM) and we explore its behavior over $\gamma_{\text{R}}$, $K$ and the particles density $\phi$. These variables were chosen specifically to investigate the thermodynamic character of the emergent structures, rather than those which derive from the strength of the self-propulsion force and external heat bath, which would together entirely determine the entropy production of a free particle without the rotational degree of freedom. For instance, MIPS is typically controlled using the Péclet number $\text{Pe}\!\propto\!v_{0}\beta\sqrt{mI\gamma\gamma_{\text{R}}}$ Bechinger et al. (2016) by varying the propulsion force, environmental temperature and relative timescales. Instead, we restrict ourselves to varying only the relative timescales through $\gamma_{\text{R}}$. Consequently, we hold all other variables constant, setting $N{=}10000$, $R{=}0.5$, $v_{0}{=}3$, $m{=}I{=}\gamma{=}1$, $\beta{=}50$ and $\epsilon{=}1$, and also utilize periodic boundary conditions.

In order to characterize the configurational change associated with MIPS we utilize the local (per particle) sixfold bond-orientational order: \big{|}q_{6}(a)\big{|}{=}\left|\frac{1}{6}\sum_{b\in{\cal N}_{a}}e^{i6\alpha_{ab}}\right|, where $\alpha_{ab}$ is the angle between $\mathbf{r}_{ab}$ and an arbitrary axis and ${\cal N}_{a}$ are the closest $6$ neighboring particles of $a$. An order parameter for the phase separation is therefore provided by the average bond-orientational order $\langle|q_{6}(a)|\rangle$. This can be complemented by statistics of the local density ${\cal X}_{d}$, defined as the empirical density within a radius $d$, since we expect a bimodal distribution under MIPS. We consider the bimodality coefficient $\zeta({\cal X}_{d}){=}(\lambda({\cal X}_{d})+1)/\kappa({\cal X}_{d})$, where $\lambda({\cal X}_{d})$ and $\kappa({\cal X}_{d})$ are, respectively, the third and and the fourth standardized moments of ${\cal X}_{d}$. The alignment within the system is instead quantified as $\rho(\theta){=}\langle 2\cos^{2}(\theta_{a}{-}\bar{\theta}){-}1\rangle$, where $\bar{\theta}$ is the mean heading across all particles. We also introduce a measure of per particle alignment $\tilde{\rho}_{a}(\mathbf{r},\theta){=}\langle 2\cos^{2}(\theta_{a}{-}\bar{\theta}_{{\cal N}_{a}}){-}1\rangle$, where $\bar{\theta}_{{\cal N}_{a}}$ is the mean heading within ${\cal N}_{a}$.

When only excluded volume interactions are considered (i.e., $K{=}0$) as expected we observe two distinct phases: a phase with MIPS and a phase without MIPS, separated by a single critical value of $\gamma_{\text{R}}$ for any given $\phi$ (see Fig. 1(a)). Analogous behavior was observed in Digregorio et al. (2018); Siebert et al. (2018). A third kinetic phase is possible when alignment interactions are included, characterized by both polar order and MIPS (see Fig. 1(b)). At density $\phi{=}0.4$, for example, this third phase is observed for values of $K\!\gtrsim\!0.006$. For lower values of $K$ the system does not exhibit polar order, however the alignment interactions affect MIPS, which occurs only at values of $K\!\lesssim\!0.002$. Importantly, the two MIPS phases with and without polar order are emergent via two distinct and incompatible mechanisms, both having distinct effects upon the thermodynamics. Explicitly, phase separation without polar order arises due to long rotational correlation times which induces jamming-like behavior, whilst phase separation with polar order arises due to flocking behavior. At intermediate $K$ there is enough alignment to reduce the correlation times of the single particle rotational dynamics, but not enough to cause global rotational correlations necessary for flocking.

The steady state, nonequilibrium, thermodynamics of the three kinetic phases is illustrated by considering two representative trajectories through the phase diagram. The first follows the onset of MIPS in the absence of alignment interactions (i.e., $K{=}0$) at fixed density $\phi{=}0.4$ by varying $\gamma_{\text{R}}$ indicated in Fig. 1(a). The structural and thermodynamic character along the trajectory is then illustrated in Fig. 2(a-b): increasing $\gamma_{\text{R}}$ up to the critical value $\sim\!0.26$ has little effect before an abrupt increase in the bimodality coefficient at the critical point indicating the onset of MIPS. This is accompanied by a decrease in mean particle velocity through jamming causing an equally abrupt change in the expected steady state entropy production. Crucially, odd and even TRS imply completely opposite variation in the entropy production rates with an even interpretation implying lower dissipation under MIPS and vice versa. This is a qualitative distinction in the thermodynamics associated with collective motion and emergent structure, not manifest in the free particle dissipation Shankar and Marchetti (2018).

The second trajectory is indicated in Fig. 1(b) for $\phi{=}0.4$ and $\gamma_{\text{R}}{=}0.3$ as MIPS without polar order is first interrupted and then reintroduced with polar order by increasing the alignment interactions through $K$. The relevant structural and thermodynamic consequences are then illustrated in Fig. 2(c-d). Polar order, measured through $\rho$, emerges beyond a critical $K{\sim}0.006$. However, spatial order is more complicated with a large and increasing bimodality coefficient abruptly dropping when the jamming mechanism is interrupted, before distinctly rising at the onset of polar order due to flocking. The bimodality coefficient then slowly increases, although not monotonically, as MIPS with polar order dominates. Below the onset of polar order the entropy production follows the spatial order as in the $K{=}0$ trajectory with mean velocity controlled by jamming. However, beyond this point the entropy production follows the polar order as the increased alignment allows for higher velocities. Once again, odd and even TRS interpretations implicate opposite variation in the nonequilibrium behavior, with the highly aligned state corresponding to high entropy production under an even interpretation.

The spatial distribution of the entropy production can be investigated by considering the dissipation associated with individual particles (see Eq. (19) and (11)). This is exemplified in Fig. 3, for the three configurations of the system previously seen in Fig. 1(g-i). In the absence of polar order the dissipative contribution from each particle closely follows the local density (Fig. 3(a), $K{=}0$ and $K{=}0.0036$). When polar order is high (Fig. 3(a), $K{=}0.1$), this trend is reversed reflecting the distinct phase separation mechanism. An odd interpretation suggests that non-polarized clusters are highly dissipative and polarized clusters are closer to equilibrium and vice versa under and even interpretation.

This ability to quantify the thermodynamic effects of specific local spatial configurations allows consideration of defects in the emergent structures. In this manner we find a nonequilibrium analogue to the increased entropies of crystalline structures due to defects. Specifically, defects are responsible for either increases or decreases in the entropy production depending on the phase and TRS interpretation. These deviations can be directly associated with individual particles. For example, Fig. 3(b) contrasts the expected entropy production rates of individual particles with their local sixfold bond-orientational order $|q_{6}|$ under MIPS without polar order. In this phase, particles along the spatial defects are characterized by higher (lower) entropy production rates compared to the particles in highly ordered regions for even (odd) $\theta$. Similarly, Fig. 3(c) contrasts the expected entropy production rates with the local alignment $\tilde{\rho}$ under MIPS with polar order. In this phase, for suitably high $K$, polar defects (as measured by $\tilde{\rho}$) are characterized by lower (higher) entropy production for even (odd) $\theta$.

In this letter we have explored the thermodynamic character that emerges from the rich collective dynamics exhibited by active matter and highlighted a hidden entropy production where rotational timescales impact dissipation in the translational degrees of freedom. Our results suggest that the richness, commonly associated with the phase structure of active matter, is mirrored in its thermodynamics, opening up a new tool to study collective phenomena on both a micro and macroscopic scale. Important questions remain, including the delicate issue of TRS which we have shown to dramatically influence any thermodynamic interpretation. We hope that the work will contribute to a deeper understanding of the thermodynamics of active systems and, more broadly, the dynamics that can lead to emergent structures.

Derivation of the entropy productions

The over-damped and under-damped models presented in the main text involve continuous Markovian dynamics described by uncorrelated stochastic differential equations (SDEs). Given the total phase space $\Omega$, the evolution of a single degree of freedom, ${x_{a}\in\{{r}_{a}^{1},{r}_{a}^{2},{v}_{a}^{1},{v}_{a}^{2},\theta_{a},\omega_{a}\}}$, $a\in\{1,\ldots,N\}$, can be expressed in terms of a deterministic component described by the function $A_{x}(\Omega(t),t)$ and a stochastic component described by the function $B_{x}(\Omega(t),t)$ in conjunction with a Wiener process:

[TABLE]

Note, here $W_{x_{a}}$ are Wiener processes satisfying $dW_{i}dW_{j}=\delta_{ij}dt$. The deterministic dynamics $A_{x_{a}}(\Omega(t),t)$ can be further divided into reversible and irreversible components Risken (1996):

[TABLE]

Diffusion coefficients can then be associated to each coordinate such that $D_{x_{a}}(\Omega(t),t)=B_{x_{a}}(\Omega(t),t)^{2}/2$.

Following Spinney and Ford (2012), noting that there is no multiplicative noise in the models, and that the dynamics in question are autonomous, we have:

[TABLE]

where the $\circ$ notation indicates a Stratonovich integration rule. We also note the convention $0/0=0$ for deterministic co-ordinates.

Entropy production in the over-damped model

First we consider the over-damped model described in the main text:

[TABLE]

Regardless of the paritity of $\theta$ under TRS, we have $A^{\text{IR}}_{\theta_{a}}{=}{\cal T}_{a}(\mathbf{r},\theta)/I$, $A^{\text{REV}}_{\theta_{a}}{=}0$, $D_{\theta_{a}}{=}(\beta I\gamma_{\text{R}})^{-1}$ and $D_{{r}_{a}^{j}}=(\beta m\gamma)^{-1}$.

Odd self-propulsion

For the odd interpretation of $\theta$ under TRS, we have $A^{\text{IR}}_{{r}_{a}^{j}}{=}-\partial_{{r}^{j}_{a}}U_{a}(\mathbf{r})/m\gamma$ and $A^{\text{REV}}_{{r}_{a}^{j}}{=}v_{0}\hat{e}^{j}(\theta_{a})$. Applying Eq. (3) we obtain:

[TABLE]

After conversion to Itō form, in the steady-state we may assume $\langle dU(\mathbf{r})\rangle{=}0$ and $\langle(\hat{e}^{j}(\theta_{a}))^{2}\rangle{=}1/2$ and thus obtain:

[TABLE]

Even self-propulsion

For the even interpretation of $\theta$, we obtain $A^{\text{IR}}_{{r}_{a}^{j}}{=}v_{0}\hat{e}^{j}(\theta_{a})-\partial_{{r}^{j}_{a}}U_{a}(\mathbf{r})/m\gamma$ and $A^{\text{REV}}_{{r}_{a}^{j}}{=}0$. Applying Eq. (3), we then obtain:

[TABLE]

Once again converting to Itō form and assuming a steady-state such that $\langle dU(\mathbf{r})\rangle{=}0$ and $\langle(\hat{e}^{j}(\theta_{a}))^{2}\rangle{=}1/2$, we obtain:

[TABLE]

Entropy production in the under-damped model

Here we consider the under-damped model described in the main text:

[TABLE]

Regardless of the odd or even interpretation of $\theta$ under TRS, we have $A^{\text{IR}}_{\omega_{a}}{=}-\gamma_{\text{R}}\omega_{a}$, $A^{\text{REV}}_{\omega_{a}}{=}{\cal T}_{a}(\mathbf{r},\theta)/I$, $A^{\text{IR}}_{\theta_{a}}{=}0$, $A^{\text{REV}}_{\theta_{a}}{=}\omega_{a}$, $A^{\text{IR}}_{{r}_{a}^{j}}{=}0$, $A^{\text{REV}}_{{r}_{a}^{j}}{=}{v}_{a}^{i}$, $D_{\omega_{a}}{=}\gamma_{\text{R}}/\beta I$, $D_{\theta_{a}}{=}0$, $D_{{v}_{a}^{j}}{=}\gamma/\beta m$, and $D_{{r}_{a}^{j}}{=}0$.

Even self-propulsion

For the even interpretation of $\theta$, we have $A^{\text{IR}}_{{v}_{a}^{j}}{=}-\gamma{v}_{a}^{i}$ and $A^{\text{REV}}_{{v}_{a}^{j}}{=}\gamma v_{0}\hat{e}^{j}(\theta_{a})-\partial_{{r}^{j}_{a}}U_{a}(\mathbf{r})/m$. Applying Eq. (3) gives

[TABLE]

Converting to Itō form, inserting the stochastic differentials $dv_{j}^{i}$ and $d\omega_{j}$ and taking expectations yields

[TABLE]

which may be straight forwardly decomposed into its per particle contributions. However, we may alternatively recognize that $\sum_{i,j}\langle v_{j}^{i}\circ dv_{j}^{i}\rangle=\langle dU(\mathbf{r})\rangle=0$ in the steady state, thus proceeding with the surviving terms to obtain

[TABLE]

allowing us to make the steady state connection

[TABLE]

Odd self-propulsion

For the odd interpretation of $\theta$, we have $A^{\text{IR}}_{{v}_{a}^{j}}{=}-\gamma{v}_{a}^{i}+\gamma v_{0}\hat{e}^{j}(\theta_{a})$ and $A^{\text{REV}}_{{v}_{a}^{j}}{=}-\partial_{{r}^{j}_{a}}U_{a}(\mathbf{r})/m$. Applying Eq. (3) we obtain:

[TABLE]

Converting to Itō form, inserting the stochastic differentials $dv_{j}^{i}$ and $d\omega_{j}$ and taking expectations yields

[TABLE]

which again forms a basis for a per particle contribution. Again, however, in the steady state we may assume $\sum_{i,j}\langle v_{j}^{i}\circ dv_{j}^{i}\rangle=\langle dU(\mathbf{r})\rangle=0$ such that we have

[TABLE]

Utilizing Eq. (17) then gives

[TABLE]

giving a starker contrast of the total contributions under odd and even interpretations of TRS for $\theta$.

Hidden entropy production beyond under-damped translational motion

The entropy production formulae for the over-damped model immediately reduce to the free particle contributions reported in Shankar and Marchetti (2018) upon setting $\partial_{r_{a}^{i}}U_{a}(\mathbf{r})=\mathcal{T}_{a}(\mathbf{r},\theta)=0$. For a single particle, the underdamped model under the same conditions for odd $\theta$ we have

[TABLE]

and for even $\theta$

[TABLE]

Thus explicit expressions depend on determining the steady state correlation $\langle\hat{e}^{i}(\theta)v^{i}\rangle$. We can calculate this explicitly in the zero inertia limit for the dynamics of $\theta$. Under such conditions we may write

[TABLE]

where $D_{\text{R}}=(\gamma_{\text{R}}I\beta)^{-1}$. By Itō’s lemma we have

[TABLE]

with analogous expression for $d\hat{e}^{2}$. This has integrating factor solution

[TABLE]

Similarly, we have an integrating factor solution for $v_{x}(t)$

[TABLE]

so

[TABLE]

or explicitly writing $\hat{e}^{1}(t_{v}^{\prime})$

[TABLE]

We compute this in the $t\to\infty$ limit corresponding to the steady state. When we do this all terms will disappear, either through the averaging, i.e.$\langle dW_{i}\rangle=0$ or through vanishing exponentials, except the term that contains $dW_{\theta}dW_{\theta}$ which we write as

[TABLE]

Sifting out with the delta correlated Wiener processes, i.e.

[TABLE]

and considering $t_{e}=t_{v}$ such that $t_{e}>t^{\prime}_{v}$ we write

[TABLE]

which becomes

[TABLE]

In the $t_{e}=t_{v}\to\infty$ limit we may safely write $\langle(\hat{e}^{1}(t^{\prime\prime}_{v}))^{2}\rangle=1/2$ and thus

[TABLE]

Computing the integral and setting $t_{e}=t_{v}=t$ we find

[TABLE]

Considering the $t\to\infty$ limit then gives

[TABLE]

and so by symmetry

[TABLE]

and thus for odd $\theta$ gives

[TABLE]

and for even $\theta$ gives

[TABLE]

in agreement with Shankar and Marchetti (2018).

The above, however, relies upon an over-damped description for the dynamics in $\theta$, not consistent with the full under-damped equation of motion. Direct computation of $\langle e^{1}(t)v^{1}(t)\rangle$ as above suffers from the non-linearity of the transform of the unit vector. However, we can utilize Eq. (17) to consider instead the long term variance $\langle(v^{1}(t)^{2})\rangle$, so as to exploit Eq. (17), which we can construct using the same integrating factor solution in Eq. (27). Expanding, taking the limit $t\to\infty$, neglecting terms that vanish in the limit and taking expectations such that first order integrals in Wiener processes vanish leaves

[TABLE]

defining $\mathcal{G}$. When paired with Eq. (17) this gives, for odd $\theta$,

[TABLE]

and for even $\theta$

[TABLE]

When $\theta$ is described by over-damped equations of motion then $\langle e^{1}(t_{1})e^{2}(t_{2})\rangle_{\text{over}}=(1/2)e^{-(I\beta\gamma_{\text{R}})^{{-}1}|t_{1}{-}t_{2}|}$ such that $\mathcal{G}=\mathcal{G}_{\text{over}}=(\gamma(\gamma+(\beta I\gamma_{\text{R}})^{-1}))^{-1}$ also in agreement with the above.

However, $\langle e^{1}(t_{1})e^{2}(t_{2})\rangle$ differs under an under-damped description. To find such a form we first consider $\langle(\theta(t_{2})-\theta(t_{1}))^{2}\rangle$. First we integrate to find $\theta(t)$

[TABLE]

from which we obtain

[TABLE]

With $\langle\omega^{2}(0)\rangle=(\beta I)^{-1}$ corresponding to the steady state we have

[TABLE]

Expecting time translation invariance and symmetry we then have

[TABLE]

Crucially, a theorem of centered Gaussian variables states William T Coffey and Yuri P Kalmykov (2012)

[TABLE]

such that we finally have (introducing subscripts for the nature of the rotational dynamics)

[TABLE]

revealing that $\mathcal{G}_{\text{under}}\geq\mathcal{G}_{\text{over}}$ and thus for odd and even $\theta$, $\langle\Delta\dot{S}^{\rm tot}\rangle_{\text{under}}\leq\langle\Delta\dot{S}^{\rm tot}\rangle_{\text{over}}$ and $\langle\Delta\dot{S}^{\rm tot}\rangle_{\text{under}}\geq\langle\Delta\dot{S}^{\rm tot}\rangle_{\text{over}}$ respectively for stationary free ABP dynamics with under-damped translational motion.

That in the steady state for a given TRS interpretation for $\theta$, coarse-graining in the rotational dynamics causes an over or under-estimation of the entropy production for any system parameters leads to the claim of a hidden entropy production associated with such coarse-graining.

The integral for $\mathcal{G}_{\text{under}}$ generally has no close form solution, but we can calculate it in the case of all free parameters set to $1$. In this case, it reads

[TABLE]

For the same parameters $\mathcal{G}_{\text{over}}=1/2$ and so the ratio $\mathcal{G}_{\text{under}}/\mathcal{G}_{\text{over}}=2e-4\simeq 1.44$, leading to the the ratios $\langle\Delta\dot{S}^{\text{tot}}\rangle_{\text{under}}/\langle\Delta\dot{S}^{\text{tot}}\rangle_{\text{over}}=2e-4\simeq 1.44$ for even $\theta$ and $\langle\Delta\dot{S}^{\text{tot}}\rangle_{\text{under}}/\langle\Delta\dot{S}^{\text{tot}}\rangle_{\text{over}}=6-2e\simeq 0.563$ for odd $\theta$.

Numerical Integration

We utilize a stochastic velocity Verlet algorithm described in detail in Kloeden and Platen (1992). For each particle, indexed by $a$, there are three momenta variables $\omega_{a}$, $v_{a}^{1}$ and $v_{2}^{2}$. Each such variable requires two zero mean, unit variance, Guassian distributed pseudo-random numbers, written $\phi_{a}^{1,1},\phi_{a}^{1,2},\phi_{a}^{2,1},\phi_{a}^{2,2},\phi_{a}^{\theta,1},\phi_{a}^{\theta,2}$, which are all mutually independent (i.e. $\langle\phi_{a}^{x,y}\phi^{m,n}_{b}\rangle=\delta_{ab}\delta_{x,m}\delta_{y,n}$). Recalling $\mathbf{r}=\{\mathbf{r}_{1},\ldots,\mathbf{r}_{N}\}$ (with similarly defined $\mathbf{v},\mathbf{\theta},\mathbf{\omega}$) and $\mathbf{r}_{a}=\{{r}_{a}^{1},{r}_{a}^{2}\}$, $\mathbf{v}_{a}=\{{v}_{a}^{1},{v}_{a}^{2}\}$, the algorithm then reads

[TABLE]

These equations were simulated using $\Delta t=0.008$.