Time-reversal symmetry violations and entropy production in field theories of polar active matter

We investigate the steady-state entropy production rate (EPR) in the Hydrodynamic Vicsek Model (HVM) and Diffusive Flocking Model (DFM). Both models display a transition from an isotropic gas to a polar liquid (flocking) phase, in addition to traveling polar clusters and microphase-separation in the miscibility gap. The phase diagram of the DFM, which may be considered an extension of the HVM, contains additional structure at low densities where we find a novel crystal phase in which a stationary hexagonal lattice of high-density ridges surround low density valleys. From an assessment of the scaling of the EPR at low noise, we uncover that the dynamics in this limit may be organised into three main classes based on the dominant contribution. Truly nonequilibrium dynamics is characterised by a divergent EPR in this limit, and sustains global time-reversal symmetry (TRS) violating currents at zero noise. On the other hand, marginally nonequilibrium and effectively equilibrium dynamics have a finite EPR in this limit, and TRS is broken only at the level of fluctuations. For the latter of these two cases, detailed balance is restored in the small noise limit and we recover effective Boltzmann statistics to lowest nontrivial order. We further demonstrate that the scaling of the EPR may change depending on the dynamical variables that are tracked when it is computed, and the protocol chosen for time-reversal. Results acquired from numerical simulations of the dynamics confirm both the asymptotic scaling relations we derive and our quantitative predictions.

The surge of interest in active matter can also be attributed to discoveries of a vast range of novel phenomena associated with motility. Particularly prominent among these are flocking and motility-induced phase-separation (MIPS). The former arises from the combined effects of activity and alignment in response to e.g. pairwise collisions in suspensions of rod-shaped particles, hydrodynamic interactions or a sensory based steering in living systems [14,19,25]. In MIPS, spherically symmetric colloidal particles interact via steric repulsion to form dense segregated clusters against a vapor background at sufficiently high average densities and long persistence times [21]. We will be primarily concerned with field theoretic formulations of dry polar flocking, where 'dry' refers to the fact that we neglect hydrodynamic interactions with the solvent fluid [19]. In addition to a local conserved density, we describe the dynamics of such flocks by a local polar order parameter. This specifies either a head-to-tail orientation of the active particles or a local swimming velocity, and breaks rotational symmetry by attaining a nonzero global value in the flocking (or polar liquid) phase. Nonetheless, we believe that many of the principles we discuss are more widely applicable -also to systems that display MIPS -and so we will view them in light of previous work that has been conducted on similar themes.
Precise identification of TRS violations from the large scale dynamics of active matter is not always trivial, as the microscopic motion does not necessarily generate global net currents. For example, in field theories of MIPS such as Active Model B (AMB), the absence of steady mass currents renders the steady-state deceivingly similar to equilibrium phase-separation [6,21]. More specifically, for phase-separating dynamics of AMB type, the local density only provides information about the underlying irreversibility through fluctuations [33]. One might ask to what extent this qualitative notion of 'looking like equilibrium' is reflected quantitatively by the entropy production rate (EPR), measuring the extent of TRS violation by the stochastic dynamics. To address this question we investigate the dominant contribution to the EPR for polar systems at low noise, allowing us to distinguish between TRS violation at the levels of fluctuating and mean global dynamics.
In this paper we observe that dynamics in the small noise regime may be organised into three main classes based on the scaling of the EPR with the strength of local noise fluctuations. Within this scheme, truly nonequilibrium dynamics is characterised by a diverging EPR in the limit of small noise. The dominant divergent contribution stems from the ground-state dynamics at zero noise, signifying the presence of steady TRSviolating currents that persist in this limit. In fact, we will show that it is possible for the dynamics to be truly nonequilibrium even when steady homogeneous mass currents do not break detailed balance alone. In this case, the violation of TRS at ground-state level is an emergent collective phenomenon which does not have any counterpart for a single active particle. When the EPR is finite in the limit of small noise, we further classify dynamics as either marginally nonequilibrium or effectively equilibrium, where the latter corresponds to the case where the EPR vanishes in this limit. Note that the EPR can also vanish on approach to a critical point while maintaining nonequilibrium behaviour but we do not address such cases here [39]. For dynamics of marginal or effectively equilibrium type, the ground-state dynamics do not violate TRS and entropy is produced only at the level of fluctuations. However, for effectively equilibrium dynamics, detailed balance is restored at small noise where we recover Boltzmann statistics to lowest nontrivial order, while it is broken by a finite amount for any infinitesimal fluctuation in the marginal case.
In the truly nonequilibrium case, the ground-state dynamics at zero noise determines the coefficient of the leading order term. When the EPR remains finite in the limit of vanishing noise, we go beyond the deterministic setting and show that we may access the leading order coefficient by including fluctuations via a systematic expansion of the dynamics in the noise strength about the steady profile. Furthermore, we confirm our predictions by explicitly comparing them with results from numerical simulations of the dynamics. We also show that the scaling exponent of the EPR can be bounded by below from symmetry arguments, and that this agrees both with the linearisation as well as simulations.
The type of field theory we consider has been granted extensive attention elsewhere in the literature. An important first analysis was performed by Toner and Tu [16,17] in their seminal approach to the hydrodynamics of flocking based on symmetry considerations of the Vicsek model. Subsequent developments included derivations of this theory via explicit coarse-graining from the Vicsek dynamics by Boltzmann-Ginzburg-Landau [9,10,40] and Dean's equation [22,28,41] approaches. For the first part of this article, we study phenomenological equations akin to those proposed by Solon et al. in [42], albeit with noise. To compare forward and time-reversed paths of this system, the local polar density must respect a discrete polar symmetry on timereversal. By following Marchetti [19] and Dadhichi [35], who consider a more general constitutive equation for the current that advects the density, we free the polar density from this constraint.
We structure the article as follows. In section 2 we introduce the Hydrodynamic Vicsek Model (HVM), where the density current is locally proportional to the polar density. Our discussion is meant to highlight the familiar phase behaviour of the model, with particular emphasis on the location of phase boundaries with respect to the phenomenological parameters of the model (as opposed to any underlying microscopic parameters). Next, in section 3 we define the EPR via the difference between timeforward and reversed path probability weights and discuss the implications of the polar density changing sign on time-reversal. Here, we also introduce the asymptotic scaling relation for the EPR at small noise, and proceed to study its behaviour in the various phases of the model. Section 4 introduces the model we refer to as the Diffusive Flocking Model (DFM), in which we consider a more general constitutive equation for the advective current that includes noise and depends nonlinearly on both density and the local polar order parameter. This allows us to consider both the case of a timeeven and time-odd polar density, and we explore how this choice changes our results from section 3. In addition, we demonstrate that in order to fully account for the entropy produced due to density currents in the flocking phase, we must also track this advective current. Finally, in section 5 we summarise our findings and present our concluding remarks and perspectives.

The Hydrodynamic Vicsek Model
Initially, the seminal numerical study by Vicsek et al. [23] inspired a large body of research on the transition to collective 'flock' motion in active particle systems with aligning interactions. Their original article considers a discrete time lattice free automaton in d = 2, where ferromagnetic spins travel in space at constant motility and align with their nearest neighbours. In essence, the transition to collective flock motion in the Vicsek model occurs due to the coupling between the XY-type spin interaction and a time-dependent connectivity matrix of spins, separating it from the classical Heisenberg model where true long range order cannot occur in d = 2 [16,17,43]. Since its inception, the model has been generalized and recast in various different forms; in continuoustime [22,28], for spatial dimensions d = 2 [12,44], to topological (rather than metric) interactions [45,46], to systems with nematic symmetry [10,12,19,47] as well as to include additional interactions such as hard-core central forces and density-dependent bare self-propulsion speeds [19,22,28].
Of the many hydrodynamic theories that have been derived from the various microscopic models, most bear resemblance with that considered initially by Toner and Tu [16,17] on phenomenological grounds, although important insights have been gained from explicit coarse-graining. For our present perspective, particularly important are those that relate the dependence of the coupling parameters in Toner and Tu's theory on the density and local polar order with the nonlinear dynamics. In particular, we adopt equations for the local particle and polarisation densities ρ and P (respectively) that are motivated by explicit coarse-graining and include the familiar microphase-separated and polar cluster regimes of the microscopic Vicsek dynamics [42,48]. Nonetheless, the approach we choose is a phenomenological one without specific reference to any underlying microscopic model.
In the following we consider a conserved density ρ(x, t) of particles, where x ∈ V and V ⊂ R 2 is a periodic domain in d = 2, which changes locally due to the flow induced by the current J and thus obeys a continuity equation Further to this, we assume that locally particles tend to move in a direction specified by the polarisation density P (x, t) = ρ(x, t)W (x, t), where W is an order parameter for the local polar order. It is important to note that we associate with P either a local head-to-tail orientation of the particles or a direction of the bare self-propulsive force.
In general therefore, one might expect the constitutive equation for the current in (1) to be some complicated expression J ≡ J (ρ, P , ∇ρ, ∇P , . . .) of ρ and P in addition to their spatial derivatives. Indeed, this will be the case for example if particles interact via steric repulsion, or if they undergo thermal Brownian motion in addition to selfpropulsion due to interactions with an underlying substrate. For now we delay this issue, which we will revisit in section 4, and focus instead on the simplest case where the constitutive equation is given by This is consistent with the case where W is the local average direction of the velocity of particles and w is the constant self-propulsion speed. Symmetry considerations alone generally lead to a high-dimensional parameter space spanned by the coefficients appearing in the equation for the polarisation density, even when the hydrodynamic expansion is truncated at lowest non-trivial order [16,17]. One of the achievements of explicit coarse-graining has therefore been to provide a more pragmatic approach to reducing the number of independent parameters by relating them to microscopic quantities. A particularly clever observation is that many of these computations lead to an equation of the form [19,35] Expressed in this way, the equation is reminiscent of a vectorial Model A (in the Halperin-Hohenberg classification [49]), with a self-advection piece, i.e. the λ-term, that explicitly violates TRS as it cannot be written as a functional derivative [33]. In fact, TRS-violation in this model is a slightly more involved issue due to the coupling between ρ and P , and we will return to this in section 3. We further take the functional F [ρ, P ] to be given by Here, F contains a local free-energy density f (ρ, P ) which is of standard quartic form, i.e.
with the notable exception that the coefficient of the quadratic term controlling the transition to the low-temperature phase depends explicitly on the local density ρ. As mentioned above, this dependence is a product of coarse-graining and is kept here in order to capture the inhomogeneous phases, while we assume all other parameters to be constant. In fact, many such procedures lead to other coefficients also carrying a nontrivial dependence on ρ, although this is the simplest dependence needed to make the isotropic-to-flock transition similar to a first order vapor-to-liquid transition. In this form, f attains the characteristic bistable form when ρ > ρ c which marks the transition to the ordered phase. In addition, in (4) we use the function Φ defined by which is often referred to as an effective pressure [9,19,40,50]. The first term on the right-hand side is the ideal gas pressure contribution. Equation (6) also implies that pressure may be reduced locally by increasing the polar order, an effect often associated with a tendency to splay in polar liquid type systems [50]. In fact, this competing effect, in which P wants to align against gradients in ρ and towards increasing |P |, culminates in an instability at sufficiently large κ leading to the formation of localized traveling polar clusters. Finally, fluctuations are accounted for in (3) via the mean-zero spatiotemporal Gaussian white noise field η with covariance where the noise coefficient D parameterises the strength of fluctuations. Importantly, the noise term is added completely ad-hoc, and is not a direct result of coarse-graining. Several authors have addressed the effects of different noise statistics, including scalar versus vectorial noise, as well as additive versus multiplicative [19,42], although all find that the phase diagram is reasonably stable against such modifications. For our purposes (7) will therefore suffice. Even in this simplified model, henceforth referred to as the Hydrodynamic Vicsek Model (HVM), we are still left with a rather large parameter space. Some reduction of this can be made by choosing suitable units; indeed we observe that under a rescaling of the time, space and the fields we may set a, b, ν and ρ c all equal to one, which we adopt in the following. Again motivated by explicit coarse-graining from microscopic Vicsek dynamics, we will also mostly be concerned with the special case in which κ = λ and w 1 = w/2, leaving us with two independent parameters (λ, w) in addition to the noise coefficient D, system size L (taking V = [−L, L] 2 for simplicity) and the conserved average local density ρ 0 = V −1 V dx ρ. With these definitions, writing out explicitly the equation for P in (3) we obtain Figure 1: Plots in a) -b) and d) -g) illustrate the ordered phases of the HVM in (1)- (2) and (8): a) -b) show two microphase-separated profiles of the system, where in a) the traveling bands form a smectic arrangement and in b) a single solitonic band travels against an isotropic background. In both plots, ρ ⊥ denotes an instantaneous (in time) average over the direction perpendicular to the motion indicated by red arrows (→). c) Phase diagram of the model at fixed ρ 0 = 1.28, κ = λ and w 1 = w/2, with data points (• , ×, , ) corresponding to figures a) -b) and d) -g) (in order of increasing λ). Solid (--) and dash-dotted (-· -) lines correspond to the phase boundaries w = λ and ρ (λ, w) = ρ 0 respectively, determined from the linear stability analysis. Figures a) -b) and d) -e) display microphase-separation (MPS), where the number of bands is seen to increase as λ is decreased. In figure f) the system is homogeneously polarized, while in g) where λ > w, both the MPS and homogeneous polar liquid (PL) phases are unstable and localized polar clusters (PC) form.
It is well known that the HVM in (1)-(3) displays both an isotropic and a polar liquid (or flocking) phase, in addition to a microphase-separated regime with both solitonic and smectic arrangements of polar bands traveling against an isotropic background [19]. As mentioned, we also find a regime in which the polar liquid state becomes unstable and traveling polar clusters emerge. In figure 1 we show typical realisations of the nonlinear steady-states and polar liquid in the HVM.
A standard linear stability analysis provides us with some insight into the nature of the phase diagram of the model, including the nonlinear phases, although we don't expect to fully classify it by such means. We provide the details of this calculation in Appendix A for completeness. In essence, one finds that the only constant and homogeneous solutions ρ 0 , P 0 to (1)-(3) at zero noise are the mean field isotropic and polarly ordered solutions. More specifically, when ρ 0 < 1 the only solution is the isotropic one for which |P 0 | = 0. When ρ 0 > 1, one additionally finds polarly ordered solutions with The isotropic phase is linearly stable only when ρ 0 < 1 beyond which global polar order emerges. However, in the miscibility gap ρ 0 ∈ (1, ρ ) where ρ (λ, w, w 1 ) = 1 + 1 2 the homogeneous polar liquid state is linearly unstable to perturbations due to the coupling between fluctuations of ρ and of P for fluctuations that are parallel to the direction of broken symmetry. In this region we observe both spatially inhomogeneous phases reported above, i.e. microphase-separation and traveling polar clusters, separated by a phase boundary that appears at sufficiently large κ. Insight into this is again provided by the linear stability analysis of the polar liquid phase, which requires that for stability. From simulations we observe that polar clusters are formed when κ > 2w 1 , both within the miscibility gap and for larger ρ 0 . We also find numerically that within the region where microphase-separation occurs, the number of bands increases with decreasing λ. Finally, note that when κ = λ, w 1 = w/2, and the inequality (11) holds, we have that ρ ∈ (5/4, 3/2). In particular, when ρ 0 is within this region, all three polarly ordered phases can be realized by varying λ and w. This is illustrated in figure 1, where we plot the resulting phase diagram for fixed ρ 0 ∈ (5/4, 3/2) in the (λ, w)-plane. All simulations were performed using a Fourier-Galerkin pseudospectral scheme with semiimplicit time stepping [51,52], initiated with both a homogeneous isotropic and polar liquid state and allowed to relax to the steady-state at several choices for D to ensure stability. There are still many aspects of the phenomenology of the HVM that deserve deeper investigation, including the exact nature of the various transitions. However, for our purposes the available knowledge is sufficient, as our analysis will only treat the stable regimes of the isotropic, polar liquid, microphase-separated and polar cluster phases. Our aim in the following will first be to investigate the EPR in the HVM, and in particular its scaling in the limit D → 0, the physical significance of which will become more readily apparent in the next section.

Entropy production at the fluctuating hydrodynamic level
Although a large body of research in statistical physics has been devoted to the study of nonequilibrium systems, arguably few general principles have emerged. Among the more important successes are the discoveries of fluctuation theorems [5]. Within the framework of stochastic thermodynamics, these have formalised the connection between entropy production and time-reversal at the level of fluctuating trajectories [4]. Informally, fluctuation theorems capture the idea that the EPR is a measure of the probabilistic disparity between observing a time-forward trajectory (or history) of a system and its time-reversal under the same ensemble. Because of this, fluctuations are essential in order to allow the time-reversed trajectory to be realisable under the time-forward dynamics. Despite this, a well-defined limit of vanishing noise strength can often be established [33,35].
Previous studies have investigated this limit of vanishing noise in field theories of active matter, e.g. for AMB describing MIPS on the hydrodynamic scale [33]. Here it was found that the scaling of the steady-state EPR at small noise depends on the phase of the system. For an isotropic system, the EPR in AMB is O(D), while it is O(D 0 ) when phase-separation has occurred. On the other hand, Dadhichi et al. noted in [35] that in their model of flocking the EPR scales as O(D 0 ) in both the homogeneous isotropic and polar liquid phases. Here we aim to provide some further details on the physics behind these results, and to organise them within a few unifying principles.
We construct the steady-state EPR from the Freidlin-Wentzell probability measure of trajectories on the time interval [−τ, τ ] [53]. For the system in (1)-(3), this is defined via the action A, where and A = ∞ otherwise. The transition probability measure P[ρ, P ] of a trajectory (ρ(t), P (t)) t∈[−τ,τ ] is then constructed in the standard way, i.e. by setting It is important to note that in this formulation of the stochastic dynamics, equation (1) acts as a constraint which limits the space of observable trajectories (i.e. those with P > 0). In order for all observable trajectories to have an observable time-reversal under P, it is necessary therefore that the protocol T we choose for time-reversal involves a polarity flip, i.e.
Indeed, this ensures that A < ∞ if and only if the composition A • T < ∞, which can be seen directly from (12). We may thus define a time-conjugate ensemble to (13) by setting which is supported on the same constrained space of trajectories as P.
Interestingly, one observes that the functional F [ρ, P ] defined in (4) is not invariant under T . In particular, under this protocol F decomposes into T -symmetric and Tantisymmetric contributions F S and F A respectively, where is the part of F that is odd in P , and F = F S + F A . In fact, because of this we cannot interpret F as a true free-energy since it would clearly have to remain invariant under time-reversal. Moreover, this also implies that the system in (1)-(3) is out of equilibrium even when λ = 0, meaning that the self-advective contribution is not the only explicitly TRS violating component in the equations of motion. Following standard treatments of stochastic thermodynamics, we formally define the steady-state entropy production rateṠ as the (log) ratio between the forward and time-reversed ensembles P and P [4, 5, 33] We assume that (17) holds in the almost sure sense. That is, we assume thatṠ = Ṡ , where · denotes a steady-state expectation, for almost all realisations of the noise under the distribution P. This assumption of ergodicity implies that we may replace noise averages by temporal averages and vice versa when computingṠ. Furthermore, definition (17) allows us to consider the EPR pathwise, i.e. as a functional of a trajectorẏ S ≡Ṡ[ρ, P ]. By construction, this functional satisfies the symmetryṠ • T = −Ṡ as can be readily observed from (17), a fact closely related with the much celebrated Gallavotti-Cohen symmetry [4]. In Appendix C we show from (12), (13) and (15) thatṠ can be expressed in integral form asṠ We will viewṠ ≡Ṡ(D) as a function of the noise coefficient D and look to determine the asymptotic scalinġ In Appendix B we show that χ ∈ {−1, 0, 1, . . .} can only take integer values. As we will argue, this result also makes sense physically. We will see that when χ = −1 the steady ground-state dynamics at D = 0 violates detailed balance in a pathwise sense, meaning that macroscopic irreversible currents are inherent to the dynamics and are not solely observable at the level of fluctuations. On the other hand, when χ = 0 the system is in a sense 'marginally nonequilibrium'. In particular, the steady D = 0 dynamics does not violate detailed balance, yet it is broken by a finite amount for any infinitesimal fluctuation and is never recovered as we send D → 0. In contrast, when χ ≥ 1, the small noise limit is effectively an equilibrium regime where detailed balance is restored. Indeed, we will show that in the isotropic phase where χ = 1 an expansion of the fields in small D allows the lowest order contribution beyond the steady D = 0 solution to be mapped onto an equilibrium system of decoupled underdamped harmonic oscillators. In the following two subsections we investigate analytically as well as numerically the scaling (19) in the various phases of the model in (1)-(3). We begin by studying the homogeneous isotropic and polar liquid phases in addition to the polar cluster phase, where some analytical progress can be made at the fluctuating level. Subsequently, in section 3.2 we look at the micophase-separated and polar cluster regimes.

Constant homogeneous ground-states
For the computations we present here, we will assume that the steady-state dynamics relaxes onto a 'ground-state' trajectory where (ρ 0 , P 0 ) solves (1)-(3) at D = 0 and that this limit is unique up to possible degeneracies arising from rotational invariance. Firstly our aim will be to classify the ground-states that satisfy the pathwise equilibrium conditioṅ In particular, if both (20) and (21) hold, the dynamics must have χ > −1. Clearly, the pathwise equilibrium ground-states include those that are invariant under T in (14), meaning that they satisfy which follows from the fact thatṠ • T = −Ṡ. The constant homogeneous isotropic state with ρ 0 = const. and P 0 = 0 provides an example of such a state. On the other hand, the polar liquid state with ρ 0 > 1 and |P 0 | = √ ρ 0 − 1 clearly violates T alone. This is where rotational invariance arises as an important symmetry principle, because it implies that P (and thusṠ) must be invariant under the parity transformation which translates to the statement that a flock is equally likely to travel to the left as to the right. Now, if the ground-state trajectory (ρ 0 , P 0 ) is PT -symmetric, i.e. it satisfies then it follows that it is pathwise equilibrium. Perhaps surprisingly then, one realises that the constant homogeneous polar liquid state in fact is pathwise equilibrium since it satisfies PT . However, this is to be expected: a charged particle gyrating at constant frequency in the plane perpendicular to an imposed magnetic field is certainly in equilibrium (although here T should be replaced by CT to include charge conjugation). Interestingly, these observations also imply that if rotational symmetry is broken a priori, for example by driving the system with an external electric field, then P would no longer be P-invariant and the polar liquid state would have χ = −1. We also note that the fact that χ > −1 in both the homogeneous isotropic and polar liquid states also follows directly from (18) by evaluating the integral at constant (ρ 0 , P 0 ). In order to go beyond the D = 0 dynamics, we must take account of fluctuations. We do so by assuming that the fluctuating dynamics admit an expansion in small √ D, following [33,35], so that Furthermore, we restrict here to case where (ρ 0 , P 0 ) is constant and homogeneous. By substituting (25), (26) into the equations of motion in (1)-(3) and collecting terms, we obtain at order D 1/2 Here, F L is the quadratic functional where we have defined the local free energy f L by in addition to the linearised effective pressure Moreover, a 0 = 1−ρ 0 +|P 0 | 2 and P 0 satisfies a 0 P 0 = 0, while η 1 is a mean-zero Gaussian white noise with We may perform a similar procedure in order to obtain an expansion ofṠ from (18) in small D of the forṁ In Appendix B we show these are the only possible terms that could enter in the expansion ofṠ, i.e. that there are no terms of half-integer order in D. Also, from the asymptotic scaling relation (19) and the positivity of the EPR, we know thatṠ k = 0 for k < χ and that the leading order coefficientṠ χ > 0. By explicitly computing this expansion ofṠ to order D 0 , it is straightforward to show that and so χ > −1 in the homogeneous phases as argued for above. For the explicit calculation of (34), we refer to Appendix C for the details. From simulations we further find that in fact χ = 1 in the isotropic phase, as shown in figure 2, although we do not explicitly computeṠ 1 . In the polar liquid case, we obtain an explicit expression forṠ 0 given bẏ In (35) we use subscripts and ⊥ to denote components of P 1 and ∇ that are parallel and perpendicular to P 0 respectively. Consistently with (34), equation (35) implies thatṠ 0 ∼ P 0 = |P 0 | for P 0 1. In fact, this could have been predicted without explicitly performing the systematic expansion ofṠ in small D. Indeed, if we were to imagine expanding the integral expression forṠ in (18) using the series representations (25), (26) at P 0 = 0, we see from simple power counting that the only combinations of fields that could possibly appear within the integrand at order D 0 are of the form Now, the symmetryṠ • T = −Ṡ excludes all of ρ 2 , ρ 2 1 and |P 1 | 2 from entering, while ∇ · P 2 would just integrate to zero over V. For the final average in (36), observe that (27) implies Hence, there are in fact no nontrivial contributions that could enter in the expansion oḟ S at order D 0 when P 0 = 0, so we must have χ ≥ 1 in the isotropic phase. SinceṠ is O(D) in the isotropic phase, we in fact recover effective equilibrium in the limit D → 0. To see this, we transform the linearised equations of motion (27), (28) to Fourier space. Throughout we use the convention that the Fourier coefficientsĥ(q) of a function h(x) are given bŷ where we with slight abuse of notation denote by V = vol(V). We thus obtain the set of equations where we have defined the damping coefficient Γ(q) = a 0 +q 2 and the noise termη 1 (q, t) is mean zero, Gaussian and white with autocovariance Polar cluster ρ 0 = 1.28, λ = 1.53, w = 1.5 Using the mappinĝ and setting X = (x, y), V = (v x , v y ) we immediately see that these follow standard equilibrium Langevin equations for an underdamped particle in a harmonic potential [2], where the potential U = ww 1 q 2 |X| 2 /2. The final degrees of freedom in (39), (40) are captured by the transverse component and q ⊥ is perpendicular to q with |q ⊥ | = q. Again this follows an equilibrium Langevin equation,V In (45), (47) the noise terms ζ, ζ T are mean zero unit variance Gaussian white noises, and interestingly the effective temperature T is defined by Since the modes V , V T are independent for all q, the dependence of the effective temperature T on q does not lead to any current in phase space. At higher order in D, however, modesρ 1 (q, t) andP 1 (q, t) are coupled at different wavevectors q via the nonlinear terms in the equation for the polar density in (3). In particular, these terms couple heat baths at different temperatures T (q), driving the dynamics at the next order away from equilibrium. Linear theory also allows us to make quantitative predictions aboutṠ 0 from (35) in the polar liquid phase. Indeed, transforming this to Fourier space we obtaiṅ where the Hermitian matrixσ pl is given bẏ In addition, in (49) we have defined the vector of Fourier modes, as well as the matrix The sum in (49) runs over modes with wavenumbers smaller than the ultraviolet cutoff Λ, which is introduced since the sum is divergent with Λ → ∞. This is sometimes seen in field theories of this kind, since they are often derived based on the assumption that they are only valid down to a certain length scale. In Appendix C we show from this thaṫ S 0 (Λ) as predicted by the linear theory diverges in the ultraviolet as Λ 2 . Although the closed form expressions for the correlators entering in (49) are too algebraically involved to report explicitly, they may be calculated straightforwardly by numerical methods. By doing this, we may calculate the corresponding sum in (49) and quantitatively compare the results with measurements ofṠ from simulations. In figure 2 we plot the results obtained from simulations, which shows good agreement with the analytical results. Plot a) in figure 2 demonstrates that the EPR is O(D) in the isotropic phase, as well as the predicted O(D 0 ) scaling in the polar liquid phase. In particular, both remain finite at fixed Λ as D → 0, withṠ → 0 in the isotropic phase.

Nonlinear ground-states: Polar clusters and microphase-separation
Previous studies have investigated the nonlinear solutions to (1)-(3) at D = 0, and particularly interesting to our present context are the seminal contributions by Solon et al. [42,48] on the structure of the banded profiles. These are effectively one-dimensional traveling wave solutions that are invariant along the direction perpendicular to the motion. We thus write P 0 = (P 0 , 0) without loss of generality, and look for solutions of the form Direct substitution then allows us to deduce a set of equations forρ 0 andP 0 in terms of the variable z = x − ct given explicitly bỹ where primes denote differentiation with respect to z. Equation (56) can be mapped onto a Newton problem for a particle in a potential under the influence of a nonlinear drag, and all stable orbits in the (P 0 ,P 0 ) plane withP 0 ≥ 0 can be uniquely identified with a pair (c, ρ g ). In terms of the stochastic equations in (1)-(3), it is assumed that the noise selects the stable steady-state profile (of which there are infinitely many [42,48]). Importantly, these solutions to (54) break both T and PT -symmetry. Thus we expect the microphase-separated steady-state to have χ = −1. Indeed, using the traveling wave ansatz in (53) and (54) we deduce two expressions forṠ −1 = lim D→0 DṠ(D), that arė The first of these is most straightforwardly derived from (C.3) in Appendix C by using the ODE (56) forP 0 , and verifies thatṠ −1 ≥ 0 as must be the case. The latter, i.e. (58), follows immediately from (18) after integrating out total derivatives. Now, from the final expression in (58) one observes thatṠ −1 vanishes identically for even distributions, i.e. those that satisfyP 0 (z 0 + z) =P 0 (z 0 − z) for some z 0 , which is exactly the PT -symmetry in (24). However, the traveling wave profiles are clearly asymmetric with a steeper front than tail, leading in general to the observedṠ −1 = 0. In figure 3 we illustrate this, and in particular how the banded profile transforms under PT . Finally, a sanity check also verifies that both expressions (57) and (58) are invariant under P alone (P 0 (z) → −P 0 (−z)) as they should be sinceṠ is, as remarked previously, oblivious to whether the wave is moving left or right (in (56) this must be complimented by c → −c). Interestingly, this mode of TRS violation at D = 0 is a collective emergent phenomenon and does not have any counterpart for a single active particle. On the microscopic scale, it is associated with different rates of promotion and demotion of spins at the head and tail of the nonlinear profile respectively, leading to a difference in the rate of dissipation from the active alignment at these edges -a point which we will explore further in a separate paper. In the following section, this point will become more explicit when we see how TRS can also be violated at D = 0 by explicitly tracking mass currents in addition to the local polar density. In particular, in this case there is an analogous mode of TRS violation on the level of a single active particle.
We include in figure 2 the scaling of the EPRṠ in both the microphase-separated and polar cluster regimes. As shown, both are truly nonequilibrium within our classification scheme, withṠ ∼ D −1 . Although we do not possess explicit polar cluster solutions to the dynamics at D = 0, it is straightforward to argue heuristically that this is what one should expect due to the highly inhomogeneous nature of the phase. For future work, we aim to investigate this in more detail.

Summary
So far we have seen that the EPR of the HVM at small noise satisfies the asymptotic By performing a small noise expansion, we may systematically investigate the coefficients that appear at each order to determine the lowest order nontrivial contribution, and thus χ. We also find that symmetries effectively bound χ from below; when the ground-state dynamics is pathwise equilibrium, the contributionṠ −1 /D =Ṡ[ρ 0 , P 0 ] is locked out and χ > −1. In the isotropic case, the fact that ρ 1 ∇ · P 1 = 0 further constrains χ > 0, and the small noise limit becomes an effective equilibrium regime. In section 4 we look at the entropy production in a generalised model, where the constitutive equation for the density current in (2) is modified to include a diffusive contribution.

TRS violations in the generalised Diffusive Flocking Model
Above we found that for the HVM, pathwise violation of detailed balance at groundstate level is the direct result of a PT -symmetry breaking by asymmetric steady D = 0 profiles, and that a steady current of density was not sufficient to cause this alone. Here we aim to show that by changing the model so as to allow independent density current fluctuations, and by tracking this current explicitly, this may no longer hold. In this case therefore, we find that the EPR in fact does diverge as D −1 due to the presence of circulating homogeneous currents of mass. Moreover, we recover an explicit expression for the pathwise EPR of a constant homogeneous polar state in this case.

The Diffusive Flocking Model
In the following, we add a diffusive contribution and noise to the constitutive equation for the density current J . Specifically, we consider J = J d + ξ as in [19], where the noise ξ is mean zero, Gaussian and white with covariance Furthermore, we take the deterministic part J d of the current J to be of the form where γ is a constant friction coefficient. Here, µ serves an analogous purpose with the chemical potential known from equilibrium diffusive systems. However, since TRS is broken in this model, there is a priori no reason that it should be the functional variation of a free-energy. Notwithstanding, we will for simplicity ignore this issue and assume that we may write with the same functional F as that which appears in the equation for P . We only make minor changes to F for stability purposes by modifying the local free-energy f (ρ, P ) to include a quadratic term in ρ, so that now In addition we add a square gradient contribution, giving us Observe, however, that these new terms do not change the equation for P in (8), since they both drop out when considering the functional variation of F with respect to P . With this choice, we have that Again, (60) is motivated by coarse graining, and the diffusive contribution arises for example in cases where interactions such as steric repulsion are included in the microscopic model [19]. Notably, there is a kind of paradigm shift when breaking the local linear relation J ∝ P , which implies that W = P /ρ should no longer be considered the local average direction of the velocity of particles. Physically, this reflects a situation on the microscopic scale where the bare self-propulsion may be thwarted by e.g. repulsive forces, so that mass currents may move against the local polar order. More importantly in the context of entropy production, this means that a trajectory of the system in which J and P do not point in the same direction is realisable in the forward time ensemble, since fluctuations alone can now reverse J at fixed P even if highly unlikely.
Numerical integration of the dynamics with J = J d + ξ, hereafter referred to as the Diffusive Flocking Model (DFM), allows us to investigate the resulting phase diagram as in section 2. On the other hand, achieving analytical progress to a comparable extent as with the HVM is more difficult. Notably, however, from a linear analysis we do in fact find a finite wavelength instability in the region where ρ 0 < 1, in which the coarsening dynamics develop a polar crystalline structure as illustrated in figure 4. Our analysis also provides us with the isotropic-to-crystal phase-boundary, and we find that in the case w 1 = w/2 the system is unstable to perturbations when (65) More specifically, when the inequality (65) holds there is a finite range q ∈ [q − , q + ] of modes that are unstable, where the exact expressions for q ± are provided in Appendix A. In contrast, we do not observe significant changes to the phase diagram in the region where ρ 0 > 1. We explain this behaviour by observing that when ρ 0 < 1, and the local polar order is weak, the diffusive dynamics is significant while it is overpowered by advective transport when the polar order is strong [19]. At D = 0 the state is stationary and has both ∂ t ρ = 0 and |∂ t P | = 0. In fact, from simulations we observe that the stronger condition | J d | = 0 is met, meaning that at D = 0 we expect From simulations we observe that the local polar order is directed such that it points in towards low density, as illustrated in figure 4. Equation (66) then tells us that the advective transport induced by P is compensated by a reversed 'diffusion' running up gradients in ρ.
In the following we also restrict to the case where D ρ = D/γ for simplicity [19], in which case the Freidlin-Wentzell action for the DFM takes the form and the path transition density P DF [ρ, P ] is constructed as before by setting P DF ∝ exp(−A DF /D). Note that in (67), we have defined the inverse gradient operator ∇ −1 = ∇ −2 ∇, i.e. with gauge choice |∇ × ∇ −1 h(x)| = 0 [33]. Crucially, with the added density current fluctuations, we are now free to define time-reversal without the polarity flip used in (14). Specifically, we let and observe that both compositions A DF • T ± are now well defined on the full space of trajectories. As before, this means that when we define the two time-reversed ensembles to P DF by setting P ± DF = P DF • T ± , all trajectories that are observable under P DF are also observable under P ± DF . Now, comparing the time-forward ensemble with each of these gives rise to two different definitions of the entropy production rate, given bẏ In section 4.2 we will attempt to understand how this choice of polar signature may alter the scaling of the EPR at low noise [35,54]. Analogously with our treatment in the previous section, we observe that when P is odd under time-reversal, the functional F splits into even and odd pieces F S and F A respectively. However, in our current setting this has further consequences as well, since it also implies that we should not consider µ a chemical potential like quantity either. Indeed, we see that µ splits into contributions Furthermore, the deterministic part of the current, J d , also splits into a P -like odd piece under time-reversal and a ∇µ S -like even piece. That is, we write Consequently, since F does not remain invariant under time-reversal, and therefore neither µ nor J d either, it could not feature in an equilibrium theory and violates TRS.
With these definitions, we may again deduce explicit integral expressions for the EPRsṠ ± by using (67)-(69) and the definitions of P DF , P ± DF , and we refer to Appendix C for the details. There we show thaṫ where K is a matrix operator with entries K αβ = ∇ −1 α ∇ β = ∇ −2 ∇ α ∇ β . Note that in (73) we have performed an average over noise histories in order to obtain the given expression, which explains why the symmetryṠ + • T = −Ṡ + does not seem to hold pathwise any longer. However, it is recovered when writing the expression out as in e.g. (C.14) in Appendix C. As in section 3, we proceed to analyse (73) and (74) in the low D limit both analytically and numerically. We also carry over the definitions we employed there, in particular defining the exponents χ ± via the asymptotic scaling relationṠ ± (D) ∼ D χ ± for D 1. As we will see, we find that similar considerations to those made before carry over in a straightforward manner, allowing us to predict the correct scaling in all cases.

Entropy production in the DFM
Continuing as in section 3, we look for ground-state trajectories (ρ 0 , P 0 ) that satisfy the pathwise equilibrium conditioṅ in order to determine when we should expect χ ± = −1. Again, it is clear that these include all states that are either T ± or PT ± -symmetric, and so the situation remains unchanged when choosing T − . Indeed, in this case the isotropic P 0 = 0 state satisfies both symmetries, while the polar liquid |P 0 | = √ ρ 0 − 1 state is PT − -symmetric only. However, we should also expect a similar situation when choosing the T + -protocol for time-reversal; since P 0 does not flip sign upon time-reversal, any constant trajectory satisfies T + alone. Thus we expect χ ± > −1 for both the isotropic gas and polar liquid, meaning that there is no clear distinction between the two protocols for the homogeneous phases at ground-state level.
On the other hand, the situation changes quite drastically once the density current dynamics are tracked explicitly. In particular, in doing so, we expect that TRS violation at the D = 0 level should become visible from a misalignment of the density current and polar density. To see this, we promote J to a dynamical variable and consider the Freidlin-Wentzell action at this level. That is, we define if ∂ t ρ + ∇ · J = 0 and A J DF = ∞ otherwise. Importantly, J now takes the role that wP had previously in section 3, in the sense that it must be odd on time-reversal. Again, this is due to the constraint imposed by the continuity equation which limits the space of observable trajectories under the action (76). Time-reversal is thus generalised accordingly by setting Pathwise there is now a clear distinction between the two protocols T ± J . Indeed, when (ρ 0 , J 0 , P 0 ) is a constant trajectory with both |J 0 | > 0 and |P 0 | > 0, the protocol T + J introduces a discrepancy between J 0 and P 0 that cannot be transformed away by parity. On the other hand, since both J 0 and P 0 transform the same way under T − J , the trajectory (ρ 0 , J 0 , P 0 ) remains invariant under PT − J as before. This fact is reflected by the EPR computed at the level of the action in (76). To see this, we set P J DF ∝ exp −A J DF /D and P J,± DF = P J DF • T ± J as before, and defineṠ ± J via the usual definition as in (69). In Appendix C we show that this computation results in the expressioṅ forṠ + J , which differs from (73) by omission of the factor K inside the first term, while in factṠ − J =Ṡ − remains unchanged from (74). It is instructive to investigate the difference between the two expressions (73) foṙ S + and (78) forṠ + J , and in particular to show that indeed the operator K = Id. To this end, we introduce the potential ϕ as the solution to the Poisson's equation which is unique up to an additive constant (assuming periodic boundaries on V). Defining P T ≡ P − ∇ϕ it follows by construction that where P T is solenoidal, i.e. ∇ · P T = 0. Importantly, this construction is in general not the same as the standard Helmholtz decomposition, since P T is not necessarily the curl of a vector potential. A specific example demonstrating that these are indeed different is provided via simplest case for which P = P 0 is constant and nonzero. Indeed, in this case, periodic boundaries forces ϕ = const. and P T = P 0 , and the latter cannot be written as the curl of a vector field that respects the periodic boundaries. The decomposition in (80) is, however, orthogonal in L 2 (V), meaning that Thus, observing that by definition we have KP = ∇ϕ, it follows that the difference betweenṠ + andṠ + J is simplẏ This is consistent with our observation that the polar liquid breaks T + J at ground-state level as remarked above, and in particular it follows that we have for the constant homogeneous ground-states. This mechanism by which TRS is broken at D = 0 due to J and P having different polar signatures under time-reversal has an analogue for a single active particle on the microscopic level. To see this we make the identificationsẋ t δ(x t − x) → J andn t δ(x t − x) → P , where x t is the position of the active particle andn t denotes its polar orientation and direction of self-propulsion. Similarly to the current J , the particle velocityẋ t must change sign under time-reversal. The polar orientationn t need not, on the other hand, and thus generally leads to a bare entropy production associated with the motility. Continuing as in section 3, we include fluctuations by performing a systematic expansion of the equations of motion and the EPRsṠ ± via the integral expressions in (73) and (74) in small D. Again, we start by assuming that the ground-state (ρ 0 , P 0 ) is constant and homogeneous. Thus, substituting the expansions in (25) and (26) into the continuity equation (1) with J = J d + ξ we obtain to lowest nontrivial order where now with slight abuse of notation In addition, the local free energy f L is given by while Φ L remains the same as in (31). The noise term ξ 1 is a mean zero Gaussian white noise process with covariance so that the linearised equation (84) is indeed independent of D. Note also that there is no change to the linearised equation for the polar density since this is the same for both models under consideration. Similarly, an expansion of the EPRsṠ ± in small D allows us to writė whereṠ ± k = 0 for k < χ ± , andṠ ± χ ± > 0 as before. By explicitly computing this expansion one finds that In particular, since the linearised continuity equation (84) of the DFM implies that the steady-state expectation ρ 1 ∇ · P 1 no longer vanishes identically as in (37), we cannot any longer expect that χ − > 0 in the isotropic phase. Since χ ± = 0 in the isotropic and polar liquid phases, we classify both phases as being marginally nonequilibrium for the DFM. This means that the linearised dynamics of the DFM cannot be mapped onto an equilibrium dynamics for any choice of P 0 . As in section 3, we now present a more in-depth analysis of the leading order term in the expansion (88) ofṠ ± . Beginning with the isotropic phase where |P 0 | = 0, we find after a Fourier transform thatṠ ± 0 may be written in bilinear form as in (49): Here, the sum runs over wavevectors q, and an explicit dependence inṠ ± 0 (Λ) on the ultraviolet cutoff Λ is introduced in order to study the limit in which it is taken to infinity. We denote the Fourier modes of ρ 1 and P 1 byρ 1 andP 1 respectively, and have definedû whereP L =q · P 1 andP T =q ⊥ · P 1 are respectively the longitudinal and transverse components of P 1 with respect toq, andq ⊥ is perpendicular to q and of unit length. The equal-time steady-state correlation matrix In addition, the Hermitian matricesσ ±,iso in (90) are given by and where we have definedw = w(1−w 1 q 2 /γw) as well as damping coefficients Γ = 1−ρ 0 +q 2 and Γ ρ = a ρ +ν ρ q 2 . In fact, we may explicitly compute C iso from the linearised dynamics and we refer to Appendix A for the details. By substituting the result of this calculation back into (90), we obtaiṅ Interestingly, we see from direct power counting thatṠ + 0 (Λ) converges as Λ → ∞, whilė S − 0 (Λ) ∼ log Λ. In fact, expressions (95) and (96) generalise trivially to dimensions d = 2, meaning thaṫ where we denote by Λ 0 a logarithmic divergence. We may perform an identical procedure in the polar liquid case, although we leave the details of this calculation in Appendix C to simplify the presentation. We note, however, that in the polar liquid phase the scaling ofṠ ± 0 (Λ) with the ultraviolet cutoff Λ changes. For our present case where d = 2, we find thatṠ + 0 ∼ Λ 2 whileṠ − 0 ∼ Λ 4 . Similarly to our treatment in section 3, we find good agreement between predictions and the results from simulations. In figure 5 we demonstrate this comparison for both of the homogeneous phases. Furthermore, all considerations extend straightforwardly to the level where we explicitly track the current J ; the scaling exponent χ + J of the EPṘ S + J is at this level given by Note, however, that in the isotropic phase the coefficient of the leading order term ofṠ + J changes by virtue of the equation (82). Finally, we consider inhomogeneous ground-states (ρ 0 , P 0 ) (or (ρ 0 , J 0 , P 0 ) at the level of J ), specifically the nonlinear polar cluster, microphase-separated and polar crystal states. Following the same reasoning as in section 3.2, we conclude that χ ± J = χ ± = −1 for both the banded profiles and polar clusters. Indeed, these profiles still break both T ± and PT ± as before and are therefore truly nonequilibrium at small noise. On the other hand, for the crystal state we find that χ + = 0, while χ − = −1, as shown in figure 5. We explain this by the observation that at D = 0 the ground-state Table 1: Scaling of the EPRsṠ ∼ D χ ,Ṡ ± ∼ D χ ± andṠ ± J ∼ D χ ± J with the noise parameter D when D 1 for the isotropic, polar liquid, microphase-separated, polar cluster and crystal regimes of the HVM and DFM.
solution is stationary, i.e. it has both ∂ t ρ 0 = 0 and ∂ t P 0 = 0. In particular, since (ρ 0 , P 0 ) is independent of time, it is in fact also invariant under T + . This may also be confirmed by inspection of e.g. (C.14), where it is apparent that the stationary condition implies we must have χ + > −1. Similarly, at the level of the current J we also find that χ + J = 0 for the crystal state. To see why this should be true, note that also |J 0 | = 0 for the polar crystal at ground-state level, meaning that there is no difference between (ρ 0 , J 0 , P 0 ) and its time-reversal under T + J for this phase. Coincidentally, this also shows that inhomogeneity is not necessarily sufficient alone to make the system truly non-equilibrium. On the other hand, the polar crystal state is clearly not invariant under T − or PT − (nor T − J , PT − J ), and so we conclude that χ − J = χ − = −1.

Conclusion
In this paper, we have studied the entropy production rate in two related models of dry polar flocks, namely the Hydrodynamic Vicsek Model and the Diffusive Flocking Model [19,35,42]. Our main results relate to the observation that the scaling of the EPR with the noise parameter D changes depending on the phase behaviour of the steadystate, and that the asymptotic scaling exponent takes integer values ≥ −1. Truly nonequilibrium behaviour is characterised by a divergent EPR in the limit D → 0, and is caused by ground-state dynamics that violate detailed balance pathwise. On the other hand, in the marginal and effectively equilibrium cases where the scaling exponent is ≥ 0, the ground-state dynamics is pathwise equilibrium and entropy is produced only at the level of fluctuations. In particular, when the scaling exponent is strictly positive, the dynamics at small noise can be mapped onto equilibrium dynamics. Both models studied display a transition from an isotropic gas to a polar liquid, in addition to nonlinear polar cluster and microphase-separated phases [42,48,50]. In the miscibility gap where microphase-separation occurs, high-density banded profiles that break parity travel against an isotropic background. For densities beyond the polar liquid threshold ρ in (10), the phase diagram is divided into two regions with a phase boundary that can be parameterised by the self-advection parameter λ and local swimming velocity w. For small λ and large w the system is a polar liquid, and as λ is increased or w decreased this state becomes unstable to perturbations leading to the formation of polar clusters. For the HVM we were able to explicitly locate both the banded-to-flock as well as the flock-to-cluster transition lines from a linear analysis. The phase diagram of the DFM, which may be considered an extension of the HVM, contains additional structure at low densities where we find a novel crystal phase in which a stationary hexagonal lattice of high-density ridges surround low density valleys. Numerical integration of the DFM also shows that the same qualitative behaviour is retained at high densities even though the density dynamics are modified by the addition of a diffusive fluctuating current. This is, however, to be expected since the diffusive dynamics are only significant to the large scale behaviour when the advective transport is comparatively small [19].
Generally for systems with polar symmetry such as those considered here, the EPR may be constructed in two different ways depending on how we choose to implement time-reversal at the level of fluctuating trajectories [35,54]. Specifically, we may choose whether the polar density should transform as a velocity-like odd quantity or a headto-tail-like even one under time-reversal, which changes the physics of the model. An exception to this is presented by the HVM, which is constructed in such a way that we only have one choice. Here, the continuity equation imposes a constraint on the space of observable trajectories, i.e. those that lie in the support of the transition probability density, which excludes the time-reversed trajectory of all observable trajectories when the polar density does not flip sign. On the other hand, when the density advection is driven by independent fluctuations, we may consider both time-signatures. In addition, we may promote the current to an explicit dynamical variable and thus construct an additional EPR at this level [33]. Surprisingly, for this latter construction, we find that the additional knowledge of the current changes the EPR only when the time-signature of the current differs from that of the polar density. When it does, detailed balance at ground-state level for a homogeneously polarised system is broken by a mismatch between the density current and polar density, analogously to the way in which timereversal symmetry may be broken on the microscopic scale by ABPs or AOUPs [54].
For both time-signatures and models considered, as well as when explicitly tracking the density current, we find that the entropy production rate diverges in the limit D → 0 in the microphase-separated and polar cluster regimes. We attribute this to the observation that both bands and polar clusters lead to traveling spatially asymmetrical profiles, which engenders a discrepancy between the time-forward and reversed movies that cannot be transformed away by parity. It is not sufficient that a profile is inhomogeneous alone, however, which is exemplified by the stationary crystal phase of the DFM. Indeed, in this case we find that when the polar density does not change sign on time-reversal, the dynamics are only marginally nonequilibrium. Also, interestingly the mode of TRS violation that causes the microphase-separated and polar cluster dynamics to be truly nonequilibrium at small noise has no analog on the microscopic scale for a single active particle, and should be considered an emergent collective phenomenon.
We also find that the polar liquid phase is marginally nonequilibrium, except in the case where we explicitly track the density current and the polar density is even under time-reversal, as noted above. When the polar density transforms like a velocity, the zero-point EPR associated with ground-state flocking vanishes due to the rotational symmetry of the dynamics. Interestingly, we may conclude from this that if we were to break rotational symmetry a priori, for example by introducing an external driving field, then flocking would in fact be truly nonequilibrium when the polar density is odd under time-reversal. This is not the case for the isotropic phase, however, which is at most marginally nonequilibrium in all cases.
We have also shown that for both the isotropic and polar liquid phases, a linearisation of the dynamics at small noise allows us access the leading order coefficient of the EPR in the marginally nonequilibrium case by evaluating steady-state averages within the linear theory. In principle, this procedure can be adapted to access coefficients at arbitrary order in an expansion in small D, although the algebra involved becomes exceedingly complex at higher orders. Moreover, we find that our analytical predictions agree well with simulations, confirming that the procedure is well suited to analyse the EPR at small noise. In table 1 we summarise the scaling of the EPRsṠ of the HVM in addition toṠ ± andṠ ± J of the DFM with the noise parameter D for the various phases investigated. equations by equating terms at O(D α ), α = 0, 1 2 , 1, . . .. Since the continuity equation is linear, we obtain the trivial hierarchy ∂ t ρ n = −w∇ · P n , n ≥ 0 (A.1) for the density coefficients ρ n . The equation for P requires more work, however. At O(D 0 ) we find that where a 0 = 1 − ρ 0 + |P 0 | 2 . Solving this equation for the polar density gives the isotropic and polar liquid solutions P 0 = 0 and |P 0 | 2 = ρ 0 − 1 respectively. At higher orders, we find that Here, F L is the quadratic functional defined in (29). The driving term ∆ P n at each order n ≥ 2 must be derived explicitly for each case, although it depends only on the fields ρ k , P k (and their spatial derivatives) for k < n. For n = 1, ∆ P 1 = η 1 is a mean zero Gaussian white noise process with covariance and in particular does not depend on D. At first nontrivial order, i.e. n = 2, the driving term is given explicitly by In particular, when P 0 = 0, we know that (ρ 1 , P 1 ) is in fact an equilibrium dynamics, albeit with Fourier modes driven by heat baths at different temperatures. The coupling of these via (A.5) consequently drives the next order process (ρ 2 , P 2 ) out of equilibrium. Moreover, since (A.3) is inhomogeneous and linear in ρ n , P n , the system of equations (A.1)-(A.3) may in principle be solved recursively to arbitrary order. Thus, we may think of the higher order driving terms ∆ P n in a similar vein. Despite this, our analysis here will be restricted to n ≤ 1.
The situation is quite similar for the DFM, although the continuity equation is no longer linear and the hierarchy in (A.1) changes accordingly. Specifically, we find that for the DFM where F L is now given in (85). Note that even though F L is modified slightly for the DFM, the hierarchy of equations in (A.3) is in fact unchanged since δF L /δP remains the same. Again, the driving term ∆ ρ n ≡ ∆ ρ n ({ρ k , P k , . . .} k<n ) is in general a nonlinear function of ρ k , P k and their gradients for k < n. For n = 1, ∆ ρ 1 = −∇ · ξ 1 , where ξ 1 is a mean zero Gaussian white noise process with covariance In the following, we wish to examine both the linear stability of the constant homogeneous ground-states as well as to deduce expressions for the correlators in the linearised theory. We will perform this calculation in two parts: first we look at the isotropic state with P 0 = 0, and subsequently the polar liquid with |P 0 | 2 = ρ 0 − 1. In both cases, we consider the more general DFM and observe that predictions for the HVM may be made by considering the limit γ → ∞.
Here,ξ L =q ·ξ 1 is the longitudinal component of the Fourier coefficientξ 1 , whilê η L =q ·η 1 ,η T =q ⊥ ·η 1 are the longitudinal and transverse components ofη 1 respectively. Because of the convention (38) we have chosen for Fourier transforms, noise correlations in Fourier space contain an extra factor of V −1 , i.e.
Observe that solutions to (A.8) are stable and decay at an exponential rate when the eigenvalues σ ± , σ T of the linear operator L iso have positive real parts. These are straightforwardly found from the characteristic equation of L iso , and are given by Equation (A.11) for σ ± may be unraveled by first understanding its behaviour at γ = ∞.
In this limit, we find that Since w, w 1 > 0 it follows that Re σ ± (γ → ∞) > 0 if and only if σ T = Γ > 0. Thus, at γ = ∞, the isotropic state is linearly stable when ρ 0 < 1, while it is unstable otherwise. From this, it is fairly straightforward to see that a similar condition holds for finite γ. Specifically, Re σ ± > 0 only when ρ 0 < 1. However, we see from (A.11) that in this case stability also requires that (A.14) After some algebra, one finds that this latter condition holds for all q if and only if Thus, the phase diagram of the DFM in the region where ρ 0 < 1 as predicted by the linear theory is no longer trivial. Indeed, when condition (A.15) is broken, there is a finite range of wave numbers q ∈ [q − , q + ] for which the corresponding modesρ 1 ,P L grow in time, where From simulations, we find that this instability leads to the polar crystal phase reported in section 4. We are also interested in calculating the equal-time correlation functions of the linear theory in order to make analytical predictions about the EPR. Since the dynamics is linear, the correlation matrix C iso in (92) solves the algebraic Riccati equation where we have defined the diffusion matrix D by To derive this, one simply needs to apply the chain rule to the left-hand side of whereû iso is defined in (51). It is well known that the solution to the Riccati equation (A.18) may be expressed in integral form. However, we find that for our present purposes it is less cumbersome to tackle it straight on. First we observe that the linear system in (A.8) reduces to the two-dimensional coupled dynamics of (ρ 1 ,P L ), in addition to the one-dimensional dynamics ofP T . Thus, clearly, we may treat these separately. Beginning with the former, the Riccati equation (A.18) may be re-expressed as a fourby-four linear system, specifically  To find the correlators of the linear theory we therefore simply invert the matrix R iso , and one may check that the solution is given by , (A.23) The final non-trivial component of C iso is found straighforwardly from the dynamics of P T , and is given by To recover the correlators in the linearised HVM, we simply take the limit γ → ∞ in (A.22)-(A.26). We obtain that In particular, this result verifies that ρ 1P * L = 0 as advertised in (37).

Appendix A.2. Polar liquid ground-state
The linearised dynamics about the homogenenous polar liquid phase may be treated similarly to the isotropic case considered above. Transforming (A.3) and (A.6) for n = 1 and |P 0 | 2 = ρ 0 − 1 to Fourier space we find that d dt where we have defined in addition to new damping coefficients As before,ξ L =q ·ξ is the longitudinal component of the noiseξ, while η and η ⊥ are respectively the parallel and perpendicular components of η 1 with respect to P 0 . Furthermore, all noise termsξ L ,η andη ⊥ are independent and of the same covariance as in the isotropic case. Rather than solving the full cubic polynomial characteristic equation of L pl , we will study its roots only for wave vectors that are parallel or perpendicular to P 0 , i.e. those for which either q ⊥ = 0 or q = 0 respectively. For both of these cases, we further assume that the limit γ → ∞ has been taken. In other words, we will study the roots of the two polynomial equations and aim to deduce the conditions under which σ ≡ σ (q ) and σ ⊥ ≡ σ ⊥ (q ⊥ ) have positive real parts. Starting with (A.32), it is straightforward to show that σ solves where the polynomial g is given by In order to determine the number of roots of g in the halfplane {Re σ > 0} we apply the argument principle (for an introduction to this method, see e.g. [55] is given by the change in the argument of g as we trace C R counterclockwise, i.e. 38) which is illustrated in figure A1. Thus, by determining Z R in the limit R → ∞, we may in particular deduce the number of roots of g with positive real part. Moreover, the choice of contour is quite convenient when dealing with polynomials such as (A.35), since we know that on the arc A R , g (Re iθ ) ∼ R 2 e 2iθ + O(R). In other words, To compute Z R in the limit R → ∞, we are therefore left with having to find the change in the argument of g along I R . Stability clearly requires that Z R → 2 as R → ∞, or combining This occurs if and only if the image g (I R ) wraps around the origin once, as illustrated in figure A1, or equivalently that the winding number of g (I R ) about the origin is 1.

g⊥(IR)
Re Im Figure A1: On the semi-circle arc A R , the polynomials g ∼ R 2 e 2iθ and g ⊥ ∼ R 3 e 3iθ . Hence, in the limit R → ∞, the number of zeros of g and g ⊥ in C R is fully determined by the winding number of g (I R ) and g ⊥ (I R ) (respectively) about the origin.
Firstly, from (A.42) we see that we must require Γ > 0, or the winding number of g (I R ) could only be 0 or −1 (recall that we are tracing the line segment I R = [−iR, iR] from iR to −iR since C R is traced counterclockwise). This is ensured so long as ρ 0 > 1, which we assume in the following. Furthermore, from (A.41), it follows that the quadratic Re g (iy) has two distinct real roots for all q = 0, given by Thus, we only need to require that Im g (iy + ) < 0 and Im g (iy − ) > 0. One may show by standard means that this occurs if and only if In conclusion therefore, it follows that all the roots of the characteristic equation of L pl are positive only when (A.44) is satisfied. Proceeding with (A.33), i.e. the characteristic equation of L pl for wave vectors that are perpendicular to P 0 , we apply the argument principle once more. In this case, the roots σ ⊥ solve the cubic polynomial equation where the coefficients c 2 , c 1 and c 0 are given by By considering the image of C R under g ⊥ , we see that along the semi-circle arc A R , the change in the argument of g ⊥ is 3πi + O(R −1 ). Thus, in this case we must require that the winding number of g ⊥ (I R ) is 3 2 in order to have Z R → 3. This occurs if and only if g ⊥ (I R ) wraps arounnd the origin in the way illustrated in figure A1. To uncover the conditions under which this occurs, we again decompose g ⊥ (iy) into its real and imaginary parts: Table A1: Leading order asymptotic expressions for the components C pl ij of the equaltime correlation matrix C pl in the limit where q → ∞, as well as when γ → ∞ is taken first.
Assuming that ρ 0 > 1, we have c 2 > 0 which is necessary in order for the winding number of g ⊥ (I R ) to be positive. Furthermore, from (A.49) we see that we must require c 0 > 0 so that the quadratic Re g ⊥ (iy) has two distinct real roots. This holds if and only if Similarly, from (A.50) we deduce that we must have c 1 > 0, or equivalently so that the cubic Im g ⊥ (iy) has three distinct real roots. Finally, requiring that Re g ⊥ (±i √ c 1 ) > 0 one straightforwardly verifies that we are indeed guaranteed that Z R → 3. This last condition is equivalent to having One may show that (A.53) holds if and only if either and For all simulations we perform, conditions (A.52) and (A.53) are satisfied. Because of this, we will only be concerned the stability requirement in (A.51). It is also worth highlighting once more that, despite the rather involved analysis undertaken here, we have not solved the full cubic characteristic equation of L pl and thus have not fully identified all necessary and sufficient conditions for linear stability. Lastly, we also investigate the structure of the correlators of the linearised theory about the polar liquid state. In analogy with the isotropic calculation, the matrix C pl in (52) solves an algebraic Riccati equation as in (A.18), although in this case all three modes (ρ 1 ,P ,P ⊥ ) remain coupled for general q. Finding its solution can be achieved by solving the linear system R pl C pl = 2D, (A. 56) where C pl = (C pl 11 , C pl 12 , . . . , C pl 33 ) T , D = (D 11 , D 12 , . . . , D 33 ) T . We avoid explicitly writing out the nine-by-nine matrix R pl here for sake of clarity of presentation, although it may be found fairly straightforwardly from the Riccati equation. Analytical inversion of (A.56) may be done using standard computer algebra systems that perform symbolic computations. Due to the algebraic complexity of the resulting expressions, we choose to only state the result in certain limiting cases. More specifically, we look at the asymptotic behaviour of the components C pl ij in the limit q → ∞ at fixed q /q ⊥ , as well as when γ → ∞ is taken first. From this, we deduce in Appendix C the dependence of the EPR on the ultraviolet cutoff Λ quoted in the main text. For the six independent components of C pl we summarize the results in table A1.

Appendix B. Structure of the EPR expansion at small noise
In this appendix, we aim to sketch a proof to show that we do not expect to see terms of order D n/2 in the expansion of the EPR. We will keep our notation in this appendix fairly general, to illustrate that this is indeed something we expect generically for field theories of this type. while u 0 is a ground-state solution to the equations of motion at D = 0. In (B.2), L is a linear operator that depends on the D = 0 solution, and we assume that its spectrum is strictly positive so that solutions to (B.2) are stable (although this is not in general sufficient for the expansion to be valid, see e.g. [57]). For the HVM and DFM, we take u = (ρ, P ) T , u n = (ρ n , P n ) T and ∆ n = (∆ ρ n , ∆ P n ) T , with In general, when the equations of motion contain multiplicative noise, ∆ 1 will also contain multiplicative factors of u 0 although this does not modify our argument. Crucially, we do however assume that ∆ 1 is mean-zero, Gaussian and white. The higher order terms ∆ n can in general be expressed as functionals of u k for k < n and ∆ 1 , i.e.
and in particular do not depend on u n . Moreover, each term that composes ∆ n must preserve order. For example, for the HVM and DFM the driving terms ∆ P 3 and ∆ P 4 can be written as linear combinations where c 1 , c 2 , c 3 , . . . and d 1 , d 2 , . . . are constants (potentially zero). More specifically, a term of ∆ n that contains α k factors of u k for 0 < k < n and β factors of ∆ 1 must satisfy Our goal in the following will be to demonstrate that the hierarchy (B.2) can be solved recursively, and that in particular the solution u n can be expressed as a linear combination of terms containing only u 0 and ∆ 1 . Then, using the fact that each term in ∆ n preserves order, we are able to determine whether the expectation of u n vanishes. More specifically, we will show that we may write To see this, note first that the general solution to (B.2) is given by Moreover, G is the linear integral operator with kernel G, which depends only on u 0 , and I is the three-by-three identity matrix. Since G is a linear operator, it follows that u n is a linear combination of terms that preserve order. The iterative solution to (B.2) is now readily found by first computing Next, we solve for u 2 and substitute in our solution for u 1 , thus Continuing iteratively, we find that , (B.14) . . .
where G k ≡ G • ∆ k . Clearly, each operation preserves order. Consequently, we have in fact shown that each u n may be written as a linear combination of terms composed only of factors of u 0 and ∆ 1 , all of order n. Since ∆ 1 is mean-zero Gaussian and white, only even moments of its distribution can be non-vanishing by Wick's theorem. Thus, it follows that Similarly, expanding the EPR aṡ it is fairly straightforward to see that each coefficientṠ k/2 , k ≥ −2, must be a linear combination of terms of order k + 2. In particular, it follows again by Wick's theorem thatṠ

Appendix C. Calculation of the EPR
Here we derive explicitly the expression (18) for the entropy production rateṠ of the HVM, as well as (73) and (74) forṠ ± of the DFM. In addition, we deduce the small noise expansion of the EPRs about the constant homogeneous isotropic and polar liquid ground-states quoted in the main text. For clarity, we choose to consider the HVM and the DFM separately. In particular, in contrast with our treatment above in Appendix A, the results obtained in the former model cannot in general be computed from the latter by sending γ → ∞.

Appendix C.1. Hydrodynamic Vicsek Model
Starting from the definitions (17) for the EPRṠ and (13) for the path transition probability density P of the HVM, one finds that we may equivalently expressṠ in terms of the T -antisymmetric part of the Freidlin-Wentzell action A, i.e.
where A = A • T as before. By applying T to A as given in (12), we find that the action A for the time-reversed ensemble may be expressed more explicitly by and A = ∞ otherwise. In particular, from (C.1) one straightforwardly deduces thaṫ In (C.3) -and throughout this appendix -we indicate by that the product with ∂ t P should be interpreted in the Stratonovich sense, i.e. employ mid-point discretisation in time.
To transform this into the expression given in (18), we simply have to observe that some terms in (C.3) are integrable. Specifically, we find that we may writė where we have defined I to be the functional i.e. as the part of F S which does not explicitly depend on the density ρ. Since we assume that the moments of P and its higher order spatial derivatives are finite in steady-state, it follows that ∆I/τ → 0 as τ → ∞ so that the first term in (C.4) may safely be ignored. In a similar vein, an integration by parts allows us to write the integral over the first term appearing in the integrand in (C.4) as where we have used the continuity equation ∂ t ρ = −w∇ · P . Again, after dividing by 4Dτ , the first term on the right-hand side of (C.6) goes away in the limit τ → ∞ since it grows sublinearly in τ . Thus, assuming that we may replace temporal averages by averages over noise realizations, we finally arrive at the expressioṅ reported in the main text. Next, we investigate the expansion of (C.7) about the constant homogeneous ground-states. To do this, we substitute the expansions (25) and (26) into (C.7) and collect terms at equal order in D. It is fairly straightforward to check that there are no contributions toṠ at orders D −1 or D −1/2 . At order D 0 we find after some fairly tedious algebra thaṫ S 0 = P 0 · V dx w (∇ · P 1 ) P 1 + ρ 1 ((2w 1 − κ)∇ (P 0 · P 1 ) + λ (P 0 · ∇) P 1 ) +2κ (∇ · P 1 ) ∇ 2 P 1 − |P 0 | 2 P 1 (C.8) In order to arrive at this expression we have only assumed that P 0 satisfies the zeroth order equation (A.2), and used (A.1) repeatedly. If we further assume that P 0 is a polarised solution with |P 0 | > 0, we may write P 0 = (P 0 , 0) without loss of generality, from which (35) follows immediately after integrating out total derivatives. After transforming (C.8) to Fourier space, we obtain (49) as stated in the main text. Using our results from Appendix A, we may now investigate the scaling of this expression with Λ → ∞. To determine this, we note that the scaling of the sum in this limit is determined by the corresponding integral in q-space, i.e.
whereσ pl is given by (50). In particular, direct substitution from table A1 and subsequently performing the integrals over q and θ yieldṡ S 0 /V ∼ P 2 0 κ 2 4π Λ 2 . (C.10) Appendix C.2. Diffusive flocking model In analogy with the above calculation, we may compute the EPRsṠ ± of the DFM from the T ± -antisymmetric part of the Freidlin-Wentzell action A DF , specificallẏ where we have defined A ± DF = A DF • T ± . Writing out the actions for the time-reversed ensembles explicitly, we find that and (C.13) Thus, proceeding as above we straighforwardly find thaṫ The expression forṠ − immediately reduces to (74) provided in the main text upon replacing the temporal average by an average over noise realizations and noting again that ∆F S /τ → 0 as τ → ∞. On the other hand, (C.14) requires a bit more massaging. Firstly, we observe that the equal-time expectation (∇ · P )∇ −2 ∇ ξ = 0, (C. 16) since here the Stratonovich product coincides with the corresponding Ito product (i.e. there is no spurious drift term). It follows that To treat the product with ∂ t P in (C.14) we must explicitly compute the spurious drift. We will show that, in fact, the spurious drift is a total derivative and therefore does not contribute to the EPRṠ + . To do this, we consider a finite discretisation of the process in Fourier space with |q| ≤ Λ, which is consistent with our numerical scheme. By applying standard stochastic calculus, we then obtain V dx [(P · ∇)P ] η = V |q|≤Λ |k|≤Λ ik β P β (q − k)P α (k) η α (−q) = D |q|≤Λ |k|≤Λ ik β (δ k, P β (k) + 2δ q,k P β (q − k)) = 0 (C.18) Here, the second equality follows from the transformation rule between Ito and Stratonovich processes [57], i.e.