Collective motion in large deviations of active particles
Yann-Edwin Keta, Étienne Fodor, Frédéric van Wijland, Michael E. Cates, Robert L. Jack
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p Collective motion in large deviations of active particles
Yann-Edwin Keta,
1, 2, 3 ´Etienne Fodor, Fr´ed´eric van Wijland, Michael E. Cates, and Robert L. Jack
1, 4 Department of Applied Mathematics and Theoretical Physics,University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom Universit´e Paris Diderot, Laboratoire Mati`ere et Syst`emes Complexes (MSC), UMR 7057 CNRS, F-75205 Paris, France D´epartement de Physique, ´Ecole normale sup´erieure de Lyon, 69364 Lyon Cedex 07, France Department of Chemistry, University of Cambridge,Lensfield Road, Cambridge CB2 1EW, United Kingdom
We analyse collective motion that occurs during rare (large deviation) events in systems of activeparticles, both numerically and analytically. We discuss the associated dynamical phase transi-tion to collective motion, which occurs when the active work is biased towards larger values, andis associated with alignment of particles’ orientations. A finite biasing field is needed to inducespontaneous symmetry breaking, even in large systems. Particle alignment is computed exactly fora system of two particles. For many-particle systems, we analyse the symmetry breaking by anoptimal-control representation of the biased dynamics, and we propose a fluctuating hydrodynamictheory that captures the emergence of polar order in the biased state.
I. INTRODUCTIONA. Motivation
Active matter emerged in the last decades as a novelclass of nonequilibrium soft systems where every con-stituent consumes and dissipates energy to produce aself-propelled motion [1–3]. It includes both living andsocial systems, such as swarms of bacteria [4, 5] birdflocks [6, 7] and human crowds [8, 9], as well as syn-thetic systems, such as vibrated particles [10, 11] andself-catalytic colloids in a fuel bath [12, 13]. In theseexperimental systems, the combination of self-propulsionand interaction can lead to collective behavior withoutany equilibrium equivalent. Collective motion with ori-entational order [10, 11] and the spontaneous formationof particle clusters despite the absence of attractive in-teractions [12, 13] are celebrated examples.Minimal models have been proposed to capture thesecollective effects, with a view to identifying the essen-tial ingredients of the dynamics which delineate genericclasses of active matter. The emergence of collective mo-tion (CM) is generally described by the Vicsek model interms of aligning active particles [14], and its equivalentToner-Tu model at hydrodynamic level [15, 16], whereaspurely repulsive active particles yielding a motility-induced phase separation (MIPS) are usually consideredto reproduce the behavior of isotropic self-propelled par-ticles [17–19]. To characterize the structure and dynam-ics, thermodynamic tools inspired by equilibrium havebeen proposed, such as pressure [20–22], others focusspecifically on the deviation from equilibrium, such asthe irreversibility of the dynamics [23–28] and the dissi-pation of energy [29–33].Several recent studies focused on large deviations of ac-tive matter [32–42]. They consider transient rare eventswhere the system does not behave ergodically. Suchevents are often accompanied by collective effects, andmay also lead to dynamical phase transitions, where atypical trajectories differ significantly from the typicalones [43–47]. Numerical techniques can be used to anal-yse these transient events by introducing a bias parame-ter which controls the distance from the typical dynam-ics [48, 49]. They open the door to studying the micro-scopic mechanisms leading to stabilize atypical collectivebehaviors. These techniques have already proved suc-cessful to unveil dynamical transitions in glassy dynam-ics [50–52] and high-dimensional chaotic chains [53, 54].It was recently shown [37] that some large deviationsof isotropic active particles are associated with CM. Inthis dynamical transition, long-ranged orientational or-der is stabilized, despite the absence of any microscopicinteractions that favour alignment. This stands in con-trast to the usual (and intuitive) expectation that CMemerges as a result of particle alignment. However, thatwork did not resolve the nature of the transition betweenthe isotropic and CM phases. This work develops furtherour understanding of this transition, including the mech-anism of spontaneous symmetry breaking, the locationof the phase transition, and the relationship of the CMwith the hydrodynamic dynamics of the system.
B. Summary of main results
Before describing our analysis, we summarise the mainresults. We consider active Brownian particles (ABPs)as a popular model of overdamped self-propelled parti-cles [17, 18]. For a long time interval of duration τ , we fo-cus on the time-averaged rate of the active work per par-ticle w τ , which quantifies how much the self-propulsionforces of particles translate into actual displacement. Theensemble-averaged rate is h w τ i . Full definitions are givenin Sec. II, below.Building on [37], we focus on large deviations where theactive work is enhanced. The resulting picture is sum-marised in Fig. 1, as a function of the active work w , andalso its conjugate field s . We restrict to situations wherethe steady state of the system is spatially homogeneous,so the activity of the particles is not enough to causeMIPS. We find that spontaneous symmetry breaking oc-curs for values of the active work w τ beyond a threshold w ∗ that is strictly greater than its average value h w τ i .There is a corresponding threshold for the biasing field,in that collective motion takes place for s < − s ∗ (thissign convention is chosen so that s ∗ >
0, it means thatthe transition takes place at s = − s ∗ and not s = s ∗ ).This result is supported by a finite-size scaling analy-sis. It resolves an open question from [37], as to whethersymmetry breaking might be present for all w τ > h w τ i ,in sufficiently large systems. For ABPs, an important pa-rameter is the rotational diffusion constant D r which de-termines the correlation time of the self-propulsion force.We find that s ∗ ∼ D r for small D r , which is the regimewhere the system differs strongly from a passive fluid.A consequence of this analysis is that the system is anisotropic fluid phase for − s ∗ < s <
0. Generic argu-ments [46, 47, 55] based on coupling between large de-viations and hydrodynamic modes mean that this phaseis hyperuniform [56]. (The CM phase may also be ex-pected to have a similar property but that question isnot addressed here.)To analyse the mechanism of symmetry breaking, wefirst solve exactly a system of two active run-and-tumbleparticles (RTPs), to demonstrate that particles naturallyalign during large deviation events. Turning to many-particle systems, we exploit connections between largedeviation theory and optimal control theory [46, 47, 57],and we also develop a Landau-Ginzburg theory for thesymmetry-breaking transition, which includes both theorientational order of the ABPs, and their hydrodynamicdensity fluctuations. This gives a detailed description ofthe CM phase.The structure of the of the paper is as follows: Sec.IIdescribes the models and outlines the theoretical back-ground; Sec. III presents numerical results for CM;Sec. IV analyses the two-particle case; Sec. V explores themechanism using large-deviation bounds based on con-trolled systems with orientational interactions; Sec. VIdescribes a Landau-Ginzburg theory for the CM transi-tion. Conclusions are summarised in Sec. VII and severalappendices contain additional technical information.
II. MODEL AND METHODSA. Active Brownian particles
We consider N active Brownian particles (ABPs) intwo spatial dimensions [37]. Their positions and orien-tations are r i and θ i . We define an orientation vector u ( θ i ) = (cos θ i , sin θ i ) which we sometimes abbreviatesimply as u i . The particles are self-propelled with (bare)speed v , they interact through a WCA interaction po- I ( w ) h w τ i w ∗ w ( s ) − s ∗ w ∗ h w τ i CM HU
FIG. 1. Schematic behaviour of (left) the rate function, fol-lowing Fig. 1(a) of [37] and Fig. 2, and (right) the active work w τ as a function of its conjugate field s , following Fig. 1(b)of [37] and Figs. 3, 4. Vertical dashed lined delimit the tworegimes identified: CM ≡ collectively moving state, HU ≡ isotropic hyperuniform state. For positive s (or equivalently w < h w τ i ), the system is phase-separated and arrested [37],that case is not discussed here. tential V ( r ) with range σ and strength ε . Define U = 1 T X ≤ i 15 the self-propulsion leads to motility inducedphase separation (MIPS) [19]. We note that since D ∝ D r , increasing ˜ l p changes the balance between the self-propulsion term and the repulsive (WCA) forces in (2).This means that large ˜ l p tends to make particles overlapmore – they appear to be softer.When presenting numerical results, we take σ = 1 asthe unit of length and we fix the time unit by setting also v = 1. For theoretical calculations, we retain v and σ as explicit quantities. B. Dissipation and active work We define the instantaneous dissipated power from apurely mechanical argument as the rate of work that theparticles exert on the solvent [58, 59]˙ W = X i ˙ r i ◦ D (cid:16) ˙ r i − √ D η i (cid:17) (4)where ◦ is a Stratonovich product. We have absorbed afactor of T into W , to obtain a reduced (dimensionless)work. Here and in the following, sums run over all par-ticles, unless otherwise stated. Using (2) and taking atime average, we find1 τ Z τ ˙ W ( t ) d t = N v D w τ + 1 τ [ U ( τ ) − U (0)] (5)where w τ = 1 v N τ X i Z τ u ( θ i ) ◦ d r i (6)is the (reduced) active work per particle [37]. This isa natural measure of how efficiently active forces createmotion. It is normalised such that h w τ i = 1 in the dilutelimit φ → 0, while h w τ i = 0 for a completely jammedsystem. For a steady state, the term involving U in (5)is zero on average, so the average dissipation is fully de-termined by the average of the active work.This active work w τ is also related to the entropy pro-duction rate in the full { r i , θ i } configuration space (whichconsiders self-propulsion as a quantity that is even un-der time-reversal[27, 37]). This differs in general fromthe entropy production measured in position space { r i } [23, 60, 61].Since (2) has three separate contributions, there is anatural decomposition of the active work w τ = 1 + w f,τ + w η,τ (7)where the constant term stems from the product of theself-propulsion direction with itself and w f,τ = − Dv N τ X i Z τ u ( θ i ) · ∇ i U d t, (8) w η,τ = 1 v N τ X i Z τ u ( θ i ) ◦ √ D η i d t . (9)On average h w η,τ i = 0, so h w τ i = 1 + h w f,τ i . The quan-tity w f,τ is negative on average, because collisions be-tween particles tend to involve particle orientation vec-tors being anti-parallel to the interparticle force. In the following, we consider situations where the system self-organises to reduce collisions, in which case w f,τ becomesless negative (it increases towards zero). C. Large deviations For any given N , the active work w τ satisfies a largedeviation principle in the limit of large τ [37] p ( w τ ) ≍ exp [ − τ N I ( w τ )] (10)where I ( w τ ) is a scaled rate function. We define thescaled cumulant generating function (SCGF) ψ ( s ) = lim τ →∞ N τ log h exp ( − sN τ w τ ) i , (11)related to I ( w τ ) by Legendre transformation (see Eq. 14below), and where we have introduced a biasing param-eter s . The SCGF can be obtained by solving an eigen-value problem, see Appendix A.There is a useful analogy between this dynamical largedeviation formalism and equilibrium statistical mechan-ics. We recall the central features of this analogy, seealso [43, 44, 47]. Trajectories of our 2-dimensional sys-tem are analogous to configurations of a 2+1-dimensionalsystem. Also, the biasing field s corresponds to a thermo-dynamic field conjugate to the active work. The SCGF ψ ( s ) corresponds to the free energy density, and is thussometimes referred to as the dynamical free energy. Anysingularity in this function is a signature of a dynamicalphase transition. In particular, we focus below on phasetransitions where rotational symmetry of the system isspontaneously broken.Continuing with this analogy, the average in (11) cor-responds to a partition function for Boltzmann-like dis-tribution of trajectories. Averages with respect to thisdistribution take the form hAi s = (cid:10) A e − sNτw τ (cid:11) h e − sNτw τ i (12)where A is a dynamical observable. Numerical compu-tation of such averages is challenging in general – it iscomparable to computing a thermodynamic average attemperature T by reweighting from an equilibrium sys-tem with temperature T ′ = T . To achieve this, we evolvesimultaneously a large population of copies of the systemto generate “biased ensembles” by cloning and deletingsome of these copies at regular steps in order to enforcethe dynamical effective Boltzmann distribution. Thismethod, known as a cloning algorithm [48, 62], allowsestimation of averages like (11,12) with a cost that scaleslinearly in τ , allowing direct access to the large- τ limit.We implement it following [49, 63], using a modified equa-tion of motion to evolve the clones. Details are given inAppendix B.Of particular interest is the quantity w ( s ) = lim τ →∞ h w τ i s (13) . . . w − h w τ i . . . . . . . D − / r I ( w ) ( w − h w τ i ) ˜ l p = 2˜ l p = 5˜ l p = 10 . 02 0 . 08 0 . . FIG. 2. Rate function I ( w ) computed with (14) rescaled by D / r for persistence lengths ˜ l p = 2 , , 10. The inset is a mag-nified version of the behavior for small w − h w τ i . Parametervalues : N = 50, φ = 0 . n c = 10 , t max = 10 . which obeys w ( s ) = − ψ ′ ( s ). Since ψ is convex [44], this isa decreasing function of s , we denote its inverse by s ( w ).The rate function I is related to the SCGF by Legendretransform, in particular I ( w ) = − ws ( w ) − ψ ( s ( w )) (14)which allows computation of the rate function from theoutput of a cloning simulation.Fig. 2 shows the rate function of the active work w ≥ h w τ i , for different persistence lengths ˜ l p . The ver-tical axis has been scaled by D − / r which leads to datacollapse near the minimum. The rate function is minimal(and equal to zero) at w = h w τ i and its curvature thereis related to the variance of the active work as1 I ′′ ( h w τ i ) = lim τ →∞ τ N [ h w τ i − h w τ i ] (15)Our data suggest that this variance is proportional to D − / r . As w increases from h w τ i , the rate function de-viates from a quadratic form, in particular its curvaturedecreases, showing that large fluctuations of w τ are lessunlikely than a simple Gaussian approximation wouldpredict. This is related to a dynamical phase transition,as we now explain. III. EVIDENCE FOR SYMMETRY BREAKINGA. Collective motion and symmetry breaking Ref. [37] focussed on a system whose parameters lie (as N → ∞ ) within the MIPS region, ˜ l p = 40 and φ = 0 . 65 (see Ref. [18] for the full phase diagram of the system).For s > i.e. biasing towards trajectories of low activework, those trajectories involve a coexistence of a densejammed, arrested domain with a dilute vapor, given thename phase-separated arrest (PSA). For s < i.e. bi-asing towards trajectories of high active work, collectivemotion (CM) is found with aligned propulsion directions,despite the absence of aligning interactions microscopi-cally. In this work, we consider exclusively trajectorieswith positive fluctuations of the active work ( s < l p , such that the unbiased behaviour of the system is thatof an homogeneous active fluid.The physical reason for CM when s < v ,they collide much less frequently, so w f,τ is increased. Inparticular, if the u i are random unit vectors then thereare large relative velocities between particles (because | u i − u j | is typically of order unity). On the other hand,perfectly aligned orientations lead to | u i − u j | = 0, sothe only sources of relative motion are the passive noises η i , η j . The larger the relative velocities, the more oftenthe particles collide, leading to smaller (more negative)values of w f,τ . Hence collective motion is a natural mech-anism for increasing w τ .It is also notable that w τ is closely related to the ratio v ( ρ ) /v that appears in theories of MIPS [22], and mea-sures the reduction in particle speed due to collisions.This further emphasises that larger active work corre-sponds to reduced collisions.Since the CM phase is associated with spontaneousbreaking of rotational symmetry, it is natural to identifyan order parameter, ν = 1 N X i u i . (16)The particle orientations u i evolve independently of theirpositions, so the steady state distribution of ν is simplythe distribution of the average of N random unit vectors.It is convenient to define also the time-average of themodulus of the order parameter: ν τ = 1 τ Z τ | ν ( t ) | dt . (17)For large times one has h ν τ i s = h| ν |i s .For large N , the central limit theorem means that p ( ν ) → ( N/π )e − N | ν | (in distribution) and hence h ν τ i ≃ r πN (18)which tends to zero as N → ∞ . On the other hand,symmetry-broken states have h ν τ i s = O (1) (19)as N → ∞ . − . − . − . . s . . . . . . h w τ i s ˜ l p = 2˜ l p = 5˜ l p = 10 − . − . − . . s con /D r . . . . . . h w τ i s / w f r ee ( s ) ˜ l p = 2˜ l p = 5˜ l p = 10 − . − . − . . s . . . . . . h ¯ ν τ i s ˜ l p = 2˜ l p = 5˜ l p = 10 − . − . − . . s con /D r . . . . . . h ¯ ν τ i s ˜ l p = 2˜ l p = 5˜ l p = 10 (a) (b)(c) (d) FIG. 3. (a) Biased average of the active work h w τ i s . (b) Thedata from panel (a) plotted using rescaled variables from(23,25), leading to data collapse. (c) Biased average of the po-larisation norm h ¯ ν τ i s (See Eq. 17) as a function of the rescaledbiasing parameter s con /D r (See Eq. 25). (d) The data frompanel (c) plotted using rescaled variables, showing data col-lapse. Parameter values (all panels) : N = 50, φ = 0 . n c = 10 , t max = 10 . As N → ∞ , the limiting value of h ν τ i s is zero through-out the isotropic phase, but non-zero in the CM phase.This leads to a singularity at the transition point s = s ∗ ,as expected for an order parameter. However, in finitesystems, the quantity h| ν |i s is always positive and has asmooth (analytic) dependence on the field s . To identifythe phase transition in numerical studies, we use that thethe finite-size scaling behaviour (18,19) is different in thetwo phases. B. Results for ˜ l p ≥ Recall that the persistence length ˜ l p measures thestrength of the active self-propulsion, compared to pas-sive diffusion. Fig. 3 shows results obtained by cloning for˜ l p ≥ 2, where the active propulsion is significant. Notethat since we focus throughout on s < 0, the point cor-responding to the unbiased (natural) dynamics is on theright of the graphs and the strength of the bias increasesfrom right to left. As the bias becomes more negative,both the active work and the orientational order param-eter increase slowly at first, before showing a more rapidincrease. Similar to (15), w ′ ( s ) = − lim τ →∞ τ N [ h w τ i s − h w τ i s ] (20)so a rapid change in w ( s ) corresponds to a large variancein the biased ensemble. The analogy with thermody-namic phase transitions suggests that the critical point − . − . − . . s con /D r . . . . . . h ¯ ν τ i s N = 50 − . − . − . . s con /D r . . . . . . h w τ i s / w f r ee ( s ) N = 50 (a) (b) FIG. 4. Biased averages of the active work h w i s [panel (a)]and of the polarisation norm h| ν |i s [panel (b)] as functionsof the rescaled biasing parameter s con /D r (See Eq. 25) fornumber of particles N = 10 , , , , Parameter values :˜ l p = 5, φ = 0 . n c = 10 , t max = 10 . s = − s ∗ coincides with the point where | w ′ ( s ) | is maxi-mal. Panels (b,d) of Fig. 3 show that appropriate changesof variable can be used to collapse the data for different˜ l p , details are given just below, in Section III C.Fig. 4 shows the dependence on system size. For − s ∗ 0, the order parameter decreases with N as in (18),while for s < − s ∗ it depends weakly on N as in (19).Together with the large variance (20), this justifies ouridentification of s = − s ∗ as a critical point at whichsymmetry is spontaneously broken. The dependence of w ( s ) on N is weaker, this function should be continuousat a critical point, with a singularity in its derivative at − s ∗ ; this is consistent with the data. C. Enhancement of self-propulsion in biasedensembles We now discuss the reasons for the data collapse ob-served in Fig. 3(b,d). Note that there are two non-trivialcontributions to w τ in (7), which affect the biased en-semble (12) in different ways. To separate their effects,observe that the average (12) in the biased ensemble maybe reformulated as hAi s = (cid:10) A e − sNτw f,τ (cid:11) v con s h e − sNτw f,τ i v con s (21)where w f,τ was defined in (8); the averages on the righthand side are computed for the natural dynamics of acontrolled ABP system in which the velocity v in (2) isreplaced by v con s = v (cid:18) − sDv (cid:19) . (22)This result was noted previously in [37, 42]. To highlightconnections between optimal-control theory and large de-viation theory [46, 47, 57], we refer generically to systemswith modified equations of motion as controlled systems ,see also Sec. V A. In this case, the only modification isthe change in self-propulsion velocity, from v to v con .Eq. (21) is an exact equality, even for finite τ . This isexplained in Appendix A, by considering the biased time-evolution operators corresponding to (12,21). Comparing(12,21), the equations of motion of the system have beenmodified, and the contribution w η,τ has been removedfrom the exponential biasing factor.Non-interacting ABPs have w f,τ = 0, in which casethis construction allows a full solution of the large-deviation problem. The average (21) in the biased en-semble reduces to an unbiased steady-state average forABPs where only the self-propulsion velocity is modified,leading to an active work w free ( s ) = 1 − sDv . (23)The linearity of this function means that the SCGF andthe rate function are quadratic: there are no collisionsso the only fluctuations of w τ come from the (Gaussian)noise, via (9).We now consider the effect of of interactions, in themodified system of (21), with propulsion velocity v con .The normalised work w f,τ is not a particularly naturalquantity in the modified system, because it contains anormalisation factor v from (8). It is more natural torescale w f by w free : the quantity that appears in theexponent of (21) is then sw f,τ = s con w f,τ w free (24)with s con = s (cid:18) − sDv (cid:19) . (25)The field s con is conjugate to the normalised work( w f,τ /w free ), it is the natural biasing parameter for thecontrolled system. Both s and s con have dimensions ofinverse time.The data collapse in Fig. 3(b,d) is obtained by plottingthe normalised work h w f,τ i s /w free ( s ) against its conju-gate variable s con , rescaled by D r . This shows that thebiased ensemble (21) is controlled primarily by the di-mensionless combination s con /D r , with a much weakerdependence on ˜ l p , at least in this regime where the ac-tive self-propulsion is strong. Hence the decompositionof w τ as (7) allows its large deviations to be analysed asa combination of two factors: the bias acts on w η,τ toincrease the self-propulsion; it acts on w f,τ to generatealignment and collective motion. The reasons why s con should be scaled by D r and not by some other rate (forexample v /σ ) will be discussed in later Sections. D. The case ˜ l p = 1 We now turn to cases where the self-propulsion isweaker, corresponding to smaller ˜ l p . Fig. 5 shows resultsfor several system sizes, comparing ˜ l p = 1 with ˜ l p = 5. − . − . . s con /D r . . . . . . . . h w τ i s / w f r ee ( s ) N = 1020 304050˜ l p = 15˜ l p = 15 − . − . . s con /D r . . . . . . . . . h ¯ ν τ i s N = 1020 304050˜ l p = 15˜ l p = 15 − . − . . s con /D r − . − . − . − . − . − . − . − . h w f , τ i s N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15˜ l p = 15 − . − . . s con /D r − . . . . . . . . . . h w η , τ i s N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15 N = 1020 3040 50˜ l p = 15˜ l p = 15 (a) (b)(c) (d) FIG. 5. (a) Biased average of the active work rescaled by thebiased free particle active work (See Eq. 23). (b) Biased av-erage of the polarisation norm h ¯ ν τ i s (See Eq. 17). (c) Biasedaverage of the force part of the active work. (d) Biased av-erage of the noise part of the active work. Parameter values : φ = 0 . n c = 10 , t max = 10 . The order parameter h ν i s no longer collapses perfectlybut the qualitative behaviour is the same as for larger˜ l p . In particular the CM transition is robust. However,the active work h w τ i s no longer collapses as a function of s con /D r , there are significant deviations. In other words,the dependence of (21) can no longer be captured by thesingle dimensionless parameter s con /D r , but it dependsalso on ˜ l p when that parameter is (relatively) small.Physically, this can be rationalised by considering twomechanisms for events with large w f,τ – in practice, theseare events where particles spend less time in contact. Thefirst mechanism is CM – if particles all move with fixedspeed v con in the same direction, they never collide, asdiscussed in Sec. III A. A second mechanism is isotropic– particles move in random directions, but tend to avoideach other when they get close. When ˜ l p is large thenthe self-propulsion velocity dominates particles’ relativevelocity, and the CM mechanism dominates. This is theregime where the data collapses in Fig. 3. We arguethat the isotropic mechanism is becoming important forsmaller ˜ l p , at which point the response to the bias s ac-quires a more complicated dependence on s/D r and ˜ l p .Note that the isotropic mechanism does not rely onactive self-propulsion and can be relevant in passive sys-tems, when considering large deviations of quantities like w f,τ . On the other hand, the CM mechanism is inherentto active systems. It is therefore not surprising that theisotropic mechanism becomes more important for small˜ l p , where the system is behaving more like a passivefluid. The limit ˜ l p → − − λLv − . − . . . τ p ψ R T P ( λ ) L/l =0 . L/l =0 . − − λLv . . . . . . ν R T P ( λ ) L/l =0 . ν RTPend ν RTPave ν RTPend ν RTPave (a) (b) FIG. 6. (a) Rescaled CGF ˜ ψ RTP = τ p ψ RTP from (C35) asa function of the rescaled biasing parameter ˜ λ = λlv , fordifferent values of the rescaled ring length ˜ L = L/l . (b) Po-larisation ν RTPave as defined in (31), showing particle alignmentfor ˜ λ < 0. The dashed lines show ν RTPend , which is evaluatedat the end of the trajectory instead of as an average over alltimes, see also (32,C64). for s con . − D r , but a more detailed understanding ofthe isotropic mechanism would be required, in order toestablish this.In addition, Fig. 5(d) shows the the noise contributionto the active work w η,τ . For ˜ l p = 1, this contribution to w τ is comparable to the force part w f,τ while for larger˜ l p , the force contribution w f,τ is the larger contribution.From (23), the noise part is of order − sD/v ; taking s ∼ D r then this scales as ˜ l − and is indeed small in theactive limit of large ˜ l p . The force part w f,τ is of orderunity in this limit. IV. SPONTANEOUS ALIGNMENT OF TWORTPS As a preliminary step before describing a detailed the-ory of CM, we illustrate the mechanism for collectivemotion by an analytic computation of large deviationsof the active work, for two RTPs on a one-dimensionalperiodic ring. Since the system is finite, there cannotbe any spontaneous symmetry breaking, but we do findthat the particles tend to align their orientations whenbiased to large active work. An outline of the large devi-ation calculations can be found in Appendix C. Analysisof the unbiased behaviour of the system can be found inRef. [64].The RTPs have positions r i for i = 1 , 2, with periodicboundaries, in a domain of size L . Their active self-propulsion velocities are α i v with α i = ± 1. Particle i tumbles with rate τ − , which corresponds to α i changingits sign, and we introduce l = v τ p the persistence length.The particles interact by a pair potential V , and thereis no thermal diffusion. Let the particle separation be r = | r − r | . Hence the equation of motion (betweentumbles) is ˙ r i = α i v − ∂∂r i V ( r ) . (26) The potential V is short-ranged with V ( r ) → ∞ as r → 0, and a length scale ǫ such that V ( r ) = 0 for r > ǫ . For particles in contact, there is a particulardistance r ∗ (less than ǫ ) such that V ′ ( r ∗ ) = − v , so thattwo particles with opposite orientations can have a force-balanced state with ˙ r i = 0. We focus on the limit of hardparticles such that ǫ → r ∗ → w RTP f = v ( α − α ) ∂∂r V ( r ) (27) w RTP f = lim τ →∞ τ Z τ ˙ w RTP f ( t ) d t (28)while there is no analogue of the noise term w η becausethere is no thermal noise. In the limit of hard particlesone has simply that ˙ w RTP f = − v if the particles are in aforce-balanced touching state [for example r = 0 with α = 1 = − α ] and ˙ w RTP f = 0 otherwise. We denoteby λ the biasing field conjugate to the active work, andintroduce the cumulant generating function ψ RTP ( λ ) = lim τ →∞ τ log D e − λ R τ ˙ w RTP f ( t ) d t E (29)such that h w RTP f i λ = − ∂ λ ψ RTP ( λ ) as usual. Here λ isplaying the role of s in the ABP system.This quantity can be obtained by solving an eigenprob-lem, as discussed in Appendix C. Fig. 6(a) shows that ψ RTP converges to a constant when λ → −∞ , there-fore lim λ →−∞ h w RTP f i λ = 0, indicating that collisions arecompletely suppressed in this regime.The analogue of (17) in this system is ν RTP τ = 1 τ Z τ α ( t ) α ( t )2 d t , (30)which is between 0 and 1, with an average value of 1 / α s are independent. Itsaverage value in the biased ensemble is ν RTPave ( λ ) = lim τ →∞ h ν RTP τ i λ . (31)Note that this quantity is averaged over the whole trajec-tory, and its determination requires one to consider boththe left- and right-eigenvectors of the associated eigen-problem. The computation is described in Appendix C,which also considers the quantity ν RTPend = (cid:28) α ( τ ) α ( τ )2 (cid:29) λ (32)which measures degree of alignment at the final time τ .Fig. 6(b) shows results. Starting from the zero-biasstate where all configurations are equiprobable ( ν RTP =1 / λ is reduced from zero,corresponding to large positive fluctuations of the activework. For large negative λ , the polarisation eventuallyreaches a plateau value. At fixed persistence length,the larger systems have weaker polarisation. On thecontrary, the polarisation decreases for positive λ , indi-cating that anti-aligned states are more probable thanaligned states, so particles spend more time in collision.As usual, the time-averaged measurement ν RTPave respondsmore strongly to the bias than the corresponding mea-surement ν RTPend at the final time [49].The main conclusion from this analysis is that a one-dimensional system of two self-propelled particles alreadyshows that biasing towards fewer collisions promotes thealignment of the particles’ orientations. We also describein Appendix C 3 a scaling regime that is relevant whenthe system is very large, which allows some simplificationof the resulting expressions. V. LARGE DEVIATION MECHANISM We now present an analysis of the mechanism by whichCM occurs in the system of many interacting ABPs, as inFigs. 2-5. For two RTPs, we have seen that alignment isa natural mechanism for suppressing collisions betweenparticles. The same is true for ABPs. To characterisethe CM state we compare the biased ensemble (12) withother ensembles that we define either by modifying theequations of motion of the system (via control forces), orby applying different kinds of bias to the system. A. Control forces Within large deviation theory, it is often useful to com-pare biased ensembles like (12) with ensembles wherethe dynamics of the system is modified by controlforces [46, 47, 57]. In principle, biased ensembles can bereproduced exactly by a suitable (optimally-controlled)system. However, these optimal control forces cannotusually be derived exactly, except in idealised models.As a generic controlled ABP system we take equationsof motion˙ r i = v con u i − D ∇ i ( U + U con ) + √ D η i , ˙ θ i = − D r ∂U con ∂θ i + p D r ξ i , (33)where U con is a control potential (dependent on all par-ticle positions and orientations) and v con is a parameter(independent of positions and orientations).For ABPs in the biased ensemble (12), we explain inAppendix A that the optimally-controlled system has v con = v con s U con = U opt s , (34)where v con s was defined in (22) and U opt s must be deter-mined from the solution to an eigenproblem, see (A2,A3).Note that the equations of motion for the optimally-controlled process have exactly the same random noise terms as the original process (2), this is a general featureof large deviations in this class of system [65].We consider several controlled systems where the ori-entational dynamics of ABPs is modified by long-rangedcoupling that favours particle alignment. We test numer-ically how well they capture the properties of the biasedensemble (12), and hence the true large deviation mech-anism. This is achieved by deriving bounds on the ratefunction I ( w ). The bounds would be exact equalities ifour approximations for the optimally controlled dynam-ics were exact. In fact the bounds are close but not exact– from this we conclude that our approximations captureimportant features of the large deviation mechanism, butthey also miss important parts of the physics. In contrastto this work, the bound considered in [42] is restricted to U con = 0, in which case the control forces only affect theself-propulsion velocity v con .To construct bounds, let P be the path probabilitymeasure for the ABPs, and let P con be the correspondingmeasure for the system with control forces. Averages inthe controlled system are denoted h . . . i con . Then for any(ergodic) controlled system with h w τ i con = w we have I ( w ) ≤ lim τ →∞ N τ D KL ( P con || P ) (35)where D KL ( Q || P ) is the Kullback-Leibler (KL) diver-gence between distributions P and Q (see Eq. A5). Inthe analogy between large deviation theory and thermo-dynamics, (35) corresponds to the Gibbs-Bogoliubov in-equality [47].As a general controlled system we take (33). In thiscase the KL divergence can be computed, see (A7,A8).The key point is that if the control forces are optimalthen (35) is an equality. We denote the path probabilitydistribution for this optimally-controlled system by P opt s so I ( w ( s )) = lim τ →∞ N τ D KL ( P opt s || P ) . (36)Averages with respect to P opt s also match averages in thebiased ensemble, that is hAi s ≃ hAi opt (see Eq. (12)), as τ → ∞ .As a general rule, the closer is the controlled system tothe true fluctuation mechanism, the more accurate willbe the bound (35). The intuition is that choosing a con-trolled process corresponds to proposing a mechanism forthe rare event, and this mechanism can occur with a par-ticular probability of order exp[ −D KL ( P con || P )]. Hence,smaller values of D KL correspond to mechanisms thatare exponentially more likely, and the mechanism thatminimises D KL is an accurate representation of the rareevent. B. Coupling among ABP orientations and upperbound on rate function As already discussed in [37], the spontaneous breakingof symmetry for s < g > U con g = − gND r | ν | (37)independent of particle positions. The direction of theorder parameter is determined by an angle ϕ through ν = | ν | (cos ϕ, sin ϕ ). The equation for the ABP orientationin the controlled system can then be written as˙ θ i = − g | ν | sin( θ i − ϕ ) + p D r ξ i . (38)Similar to the original ABPs, this orientational equationof motion is independent of all particle positions. Hence,integrating out the particle positions leads to a mean-fieldsystem of interacting rotors, which is fully described by(38). For large N , appendix D 1 shows that this systemspontaneously breaks rotational symmetry at g = D r .That is, for g > D r then h| ν |i con = O (1) as N → ∞ , but h| ν |i con = O ( N − / ) for g < D r . The resulting situationis similar to (19,18).It was argued in [37] that this very simple controlledmodel can already capture quite accurately the collec-tive motion phase, and it can predict the rate functionin this regime. To explore this idea in more detail, wewrite P con g for the path probability distribution for theABP dynamics with control potential U con g . Using thisdistribution with Eq. 35 yields an upper bound I ( w con g ) ≤ lim τ →∞ N τ D KL ( P con g || P ) (39)where w con g = h w τ i con with control potential U con g . Thisbound would be an equality if the controlled model fullycaptured the CM phase. It will be compared with theexact result in Sec. V D. C. Large deviations of order parameter and lowerbound on rate function We also derive a lower bound on the rate function,similar to [37]. To achieve this, we consider large devi-ations of the orientational order parameter. Specificallywe consider ν τ , defined in (17) as the time-average of themodulus of the order parameter. The statistical proper-ties of other similar quantities (for example the modulusof the time average) have different dependence on N, τ ,so some care is required in the following arguments.The quantity ν τ obeys an LDP as τ → ∞ which wewrite as p ( ν τ ) ≍ exp [ − τ N J ( ν τ )] . (40) where J is the rate function. For large N , the func-tion J can be obtained analytically, see Appendix D. Inparticular, one has for large N and small ν thatlim N →∞ J ( ν ) = 12 D r ν + O ( ν ) . (41)Moreover, in this joint limit of large N and small ν , theoptimally controlled dynamics associated with these largedeviations can be captured exactly by (38). The largedeviations of ν also satisfy a (general) bound analogousto (35), which is J ( ν ) ≤ lim τ →∞ N τ D KL ( P con || P ) (42)which holds for any controlled dynamics with h| ν |i con = ν . These results can be used to obtain a bound on I ( w ). We write ν ( s ) = h| ν |i s . Now consider (42) with P con = P opt s , as the optimally-controlled dynamics forlarge deviations of the active work, as in (36). Thisoptimally-controlled dynamics has h w τ i con = w ( s ) and h| ν |i con = ν ( s ). Combining (42) with (36), we obtain I ( w ( s )) ≥ J ( ν ( s )) . (43)This is a lower bound on I ( w ), which was derived in [37]by the contraction principle of large deviations. From(42), the bound (43) is exact if the optimally controlleddynamics for large deviations of w (that is P opt s ) is alsoan optimally-controlled dynamics for large deviations of ν . Note that evaluation of (43) requires knowledge of ν ( s ),which is not available analytically – instead one mustperform cloning simulations. For this reason, (43) is nota predictive result. However, it is a useful result becausethe accuracy of the bound reveals the extent to which the(unknown) mechanism for large deviations of w is similarto the (known) mechanism for large deviations of ν , aswe now discuss. D. Numerical evaluation of bounds Fig. 7 shows results for the rate function I ( w )(Sec. II C), compared with the upper and lowerbounds (39,43). We show data for ˜ l p = 5 as well as˜ l p = 40, which was the case considered in [37].To evaluate the upper bound, we perform unbiasedsimulations of ABPs with the rotational equation of mo-tion given by Eq. 38, over a range of torque parameters g . We compute D KL ( P con g || P ) using (A10) and the aver-age active work h w τ i con . The upper bounds in Fig. 7 areparametric plots using these data (with g as the param-eter).For the lower bound, we compute J (¯ ν ) from cloningsimulations of rotors (see Appendix D 3) and ν ( s ) fromcloning simulations of ABPs, then compose these func-tions to obtain the right hand side of Eq. 43. We stress0 . 00 0 . 05 0 . 10 0 . 15 0 . 20 0 . 25 0 . 30 0 . 35 0 . 40 0 . 45 0 . w − h w τ i . . . . . . . . . . . . . I ( w ) ˜ l p = 5 I ( w ) J ( h ¯ ν τ i s ( w ) ) D ( P con g ( w ) || P ) . . . . . . w − h w τ i − . − . − . . . . . . B ( w ) / I ( w ) − B ≡ J B ≡ D . 00 0 . 05 0 . 10 0 . 15 0 . 20 0 . 25 0 . 30 0 . w − h w τ i . . . . . . . I ( w ) ˜ l p = 40 I ( w ) J ( h ¯ ν τ |i s ( w ) ) D ( P con g ( w ) || P ) . . . . w − h w τ i − . − . − . . . . . . B ( w ) / I ( w ) − B ≡ J B ≡ D (a) (b) FIG. 7. (a) Rate function I ( w ) from the cloning algorithm at ˜ l p = 5, compared with the upper and lower bounds (39,43).For the rate function, the solid line is a spline interpolation. To evaluate bounds, results for h ¯ ν τ i s ( w ) and J (¯ ν ) are obtainedby the cloning algorithm, and D ( P con g ( w ) || P ) is obtained by simulation of the controlled dynamics. Fourth order polynomialinterpolation is then used to obtain values at the desired values of w . (inset) Relative error between the bounds and the ratefunction. (b) Similar results as (a), for ˜ l p = 40. General parameter values : N = 10, φ = 0 . Cloning parameter values : n c = 10 , t max = 10 , N runs = 10. Parameters for controlled dynamics : t max = 10 , N runs = 10. that the number of particles N and the rotational diffu-sivity D r have to be consistent between simulations forthis comparison to hold.The bounds of Fig. 7 capture the main features ofthe rate function but there are significant deviations, es-pecially for smaller ˜ l p . We note in particular that (i)the bounds (39,43) are not accurate for small values of w −h w i but become accurate for larger w ; (ii) the boundsare more accurate for larger ˜ l p . We now discuss these ob-servations. E. CM (symmetry-broken) state From Fig. 3, the system is in the collective motionphase for w > w ∗ with w ∗ ≈ . 4, weakly dependent on˜ l p . For ˜ l p = 5, the CM phase is w − h w τ i & . 15. In thisrange, the inset of Fig. 7 shows that upper and lowerbounds are accurate to around 20% of the rate function.For ˜ l p = 40, the accuracy of the lower bound is evenbetter.In both cases, the lower bound is more accurate. Partof this effect may be attributed to the fact that the up-per bound does not account for the enhancement of self-propulsion in the biased ensemble (Sec. III C). The lowerbound does account (at least partly) for this effect, sinceit uses the functions w ( s ) and ν ( s ) as computed in cloningsimulations. The upper bound could be improved by in-cluding the controlled velocity v con as an additional vari-ational parameter in (35) and optimising it numerically,similar to [42, 66]. However, it is sufficient for our argu- ment to keep v con = v .The conclusion for this regime is that fluctuations ofthe active work are strongly coupled to those of the orien-tational order parameter. As a result, the bounds (39,43)can capture the behaviour of the rate function almostquantitatively. As noted in [37], the log-probability of alarge fluctuation of w τ is almost the same as that of thecorresponding orientational fluctuation. F. Isotropic state h w τ i < w < w ∗ When the active work is close to its average value,the system does not break symmetry (recall Fig. 3) andone also sees that the bounds in Fig. 7 do not capturethe rate function in an accurate way. The symmetry-breaking transition happens at w = w ∗ . As N → ∞ ,isotropic systems (with w < w ∗ ) have ν ( s ) = 0 so thelower bound (43) tends to zero. The rate function is notzero so the bound is not at all accurate in this range. Theconclusion for this regime is that fluctuations where w τ is enhanced can also occur by an alternative fluctuationmechanism (not CM), and that this mechanism is domi-nant in the isotropic phase. As noted above, the modifiedself-propulsion v con = v con s is relevant in this regime, andparticles may also be repelled from each other withoutany alignment, as discussed in Sec. III D. We explain inSec. VI that hydrodynamic density fluctuations are alsorelevant. In other words, several effects act to enhancethe active work in this regime, there is no collective mo-tion, and the bounds of Fig. 7 do not follow the rate1 − . − . − . . s/D r . . . . . . . . h w τ i s N =10 2050˜ l p = 5˜ l p = 5 . 25 0 . 45 0 . h w τ i . . . . . . . . . h ¯ ν τ i N = 10 N = 20 N = 50 ˜ l p = 2 h . . . i s h . . . i g h . . . i s h . . . i g − . − . − . . − g/D r . . . . . . . . h w τ i g ˜ l p = 5 N =10 2050 N =10 2050 . 25 0 . 45 0 . h w τ i . . . . . . . . . h ¯ ν τ i N = 10 N = 20 N = 50 ˜ l p = 5 h . . . i s h . . . i g h . . . i s h . . . i g − . − . − . . s/D r − . . . . . . . . . . g ( w ( s )) / D r N =10 2050˜ l p = 5 − s/D r ˜ l p = 5 − s/D r . 25 0 . 45 0 . h w τ i . . . . . . . . . h ¯ ν τ i N = 10 N = 20 N = 50 ˜ l p = 10 h . . . i s h . . . i g h . . . i s h . . . i g (a) (b)(c) (d)(e) (f) FIG. 8. (a) Active work as a function of the biasing pa-rameter from cloning simulations. Solid lines are spline in-terpolation. (b) Parametric plot of h| ν |i = (cid:10) τ R τ d t | ν ( t ) | (cid:11) as a function of the active work. We show results for bi-ased ensembles ( h . . . i s ) and for controlled dynamics ( h . . . i g ). (c) Active work as a function of the torque parameter w inthe controlled (modified) dynamics. Solid lines correspondto spline interpolations, points correspond to numerical data. (d) Similar to (b), with ˜ l p = 5. (e) Composition of the re-lation g ( w ) from modified dynamics and w ( s ) from cloning. (f) Similar to (b,d), with ˜ l p = 10. General parameter val-ues : φ = 0 . Cloning parameter values : N clones = 10 , t max = 10 , N runs = 10. Controlled dynamics parameter val-ues : t max = 10 , N runs = 10. function accurately. G. Comparison with controlled system To conclude our discussion of orientational fluctua-tions, we consider Fig. 8 which further illustrates therelationship between controlled dynamics and the large-deviation mechanisms. Panels (a,c) show the dependenceof the active work, as either s or g is varied. Also, panels(b, d, f) show the orientational order parameter.Note that s is a biasing field in the sense of large devi-ations while g is physical coupling between ABP orienta- tions – nevertheless, the response to both kinds of pertur-bation is similar, in that both h w τ i λ and h| ν |i λ increase.The responses are different when the system is close toits unbiased steady state, in that the biased ensemble re-sponds mostly by a change in h w τ i λ while the controlledsystem responds mostly by increasing h| ν |i λ . However,when the system is perturbed further from its steadystate, both systems respond by entering a CM state, inwhich their behaviour is similar, with long-ranged orien-tational order and enhanced active work.Finally, Fig. 8(e) compares how large a coupling g isrequired in order to achieve an active work equal to w ( s ).In the CM phase, one sees that g ( w ( s )) ≈ − s . From(D5), rotational symmetry is spontaneously broken forthe controlled system with g > D r . Hence Fig. 8(b, d,f) are consistent with the observation from Fig. 3 thatsymmetry is broken in the biased ensemble for s < − s ∗ with s ∗ ≈ D r . A concrete theoretical explanation for thiseffect will be given in the next Section. VI. LANDAU-GINZBURG THEORY We have explained that many aspects of the collectivemotion phase of ABPs in terms of mean-field interac-tions among their orientations. In this section we discusshow this mean-field picture can be embedded in a hy-drodynamic theory, similar to Macroscopic FluctuationTheory [67]. A. Theory The minimal description that allows modelling of thesymmetry-broken state is to consider a local density field ρ and a corresponding polarisation field P . These are de-fined on hydrodynamic length scales: given a mesoscopicregion Ω r centred at r , we define ρ ( r ) = 1 | Ω r | Z Ω r X i δ ( r − r i )d r (44) P ( r ) = 1 ρ ( r ) | Ω r | Z Ω r X i u i δ ( r − r i )d r (45)where | Ω r | denotes the volume of Ω r . The polarisation P is normalised as the average orientation (so | P | < ρ is a slow hydro-dynamic field, in the sense that density fluctuations onlarge length scales ℓ relax on long time scales of order ℓ . On the other hand P is a fast field, in that polari-sation fluctuations on all scales relax quickly (on a timescale of order D − r ) to quasi-steady states which dependin general on ρ .As a minimal model for ABPs in biased ensembles,we propose (in a “top-down” coarse-grained approach)the following (It¯o) equations of motion for ( ρ, P ), similar2to [68]:˙ ρ = −∇ · JJ = J d + p σ ( ρ ) η , ˙ P = − γ ( ρ, P ) f ( P ) + b ( ρ, P ) ∇ ρ + p γ ( ρ, P ) ξ (46)where σ and γ are noise strengths; b is a coupling betweenpolarisation and density gradients; f is a thermodynamicforce that acts on the polarisation; J d is the deterministicpart of the current; and η is a Gaussian white noise withzero mean and variance h η α ( t ) η β ( t ′ ) i = δ αβ δ ( t − t ′ ). Wetake a deterministic current J d = v ρ P − D c ( ρ ) ∇ ρ (47)where the first term incorporates the effect of self-propulsion, while D c ( ρ ) is the hydrodynamic (collective)diffusion constant, which depends on density. Our the-ory is restricted to states without MIPS, this requiresthat D c ( ρ ) > ρ . (In fact the diffusion constant isadditionally renormalised by polarisation fluctuations, asdiscussed in Appendix E, and it is the renormalised dif-fusion constant that should be positive.) In the absenceof MIPS, it is consistent to assume that ∇ ρ is O ( ℓ − )where ℓ is the hydrodynamic length scale. This is thereason that higher gradients of ρ are neglected in (46).It is also useful to compare (46) with [69], which is arigorous hydrodynamic theory for a similar model, ex-cept that the polarisation is a slow field in that case. Inthe ABP context, their analysis corresponds to reduc-ing v and D r as the system size increases, so that thepolarisation becomes a hydrodynamic variable. This cor-responds to an idealised limiting case, but it shows howsuch theories can be justified rigorously. Comparing (46)with the Toner-Tu theory [15] and other theories for CMsuch as [70], our minimal model (46) has fewer couplings,in particular it lacks advective terms in the polarisationequation. Such terms are not relevant at the level con-sidered here, but might be required for a quantitativedescription of fluctuations in the CM phase.On the hydrodynamic scale, it is consistent to approxi-mate the (total) active work of a trajectory as a functionof the density and polarisation fields: w τ N τ = Z τ Z [0 ,L ] ω ( ρ, P ) d r d t (48)where ω ( ρ, P ) is the typical (average) active work perunit volume, in a region with density ρ and polarisation P .Now consider a large system of linear size L , and definehydrodynamic co-ordinates as˜ r = r /L , ˜ t = Dt/L . (49)We also express the fields in these rescaled variables as˜ ρ (˜ r , ˜ t ) = ρ (˜ r L, ˜ tL /D ), and ˜ P (˜ r , ˜ t ) = P (˜ r L, ˜ tL /D ).[Note, there is no rescaling of the fields themselves, hence | ˜ P (˜ r , ˜ t ) | ≤ N − ˜ l p =2 510 τ ( N ) =0 . N . N . Nτ ( N ) =0 . N . N . N V a r( ν τ ) − C o v ( | ν τ | , w τ ) N . . . . . . N τ ( N ) V a r( w τ ( N ) ) ˜ l p = 5 τ ( N ) = 2 . N ˜ l p = 5 τ ( N ) = 2 . N N . . . . . h w τ i (a) (b) FIG. 9. (a) Normalised covariance of the active work w τ andthe squared average polarisation ν τ = (cid:0) τ R τ ν ( t ) d t (cid:1) . (b) Normalised variance of the active work with a log N fit indashed line. Parameter values : φ = 0 . N runs = (a) · (b) · . and ˜ ρ (˜ r , ˜ t ) σ π/ φ = ¯ ρσ π/ ρ = ( ∂ ˜ ρ/∂ ˜ t ). We define ˜ J = J L/D so that thecontinuity equation retains its form:˙˜ ρ = −∇ · ˜ J . (50)Using (46) and working in hydrodynamic co-ordinates,we obtain a formula for the probability of a trajectory inthe biased ensemble,Prob (cid:2) ˜ ρ, ˜ J , ˜ P (cid:3) ∝ exp − L Z ˜ τ Z [0 , S (cid:2) ˜ ρ, ˜ J , ˜ P (cid:3) d˜ t d˜ x ! (51)which is valid only if (50) holds (else the probability iszero); and the Lagrangian is S = D σ | ˜ J − ˜ J d | + 14 Dγ (cid:12)(cid:12)(cid:12)(cid:12) DL ˙˜ P + γ f L − b ∇ ˜ ρ (cid:12)(cid:12)(cid:12)(cid:12) + sL D ω (52)where ˜ J d = J d L/D and we have omitted the dependenceof σ, γ, ω on (˜ ρ, ˜ P ), for ease of writing.The terms in (52) that involve J and ω are familiarfrom other analyses of large deviations in systems withhydrodynamic modes [55, 71]. In particular, the fact that s enters through the combination sL is familiar fromearlier studies [72], it reflects the fact that hydrodynamicdegrees of freedom respond strongly to the bias, becauseof their slow relaxation. However, the factors of L thatappear in the term involving ˜ P are not expected in hy-drodynamic theories, they reflect the fact that P is a fastfield.3 B. Mean-field analysis We first consider the case where (˜ ρ, ˜ J , ˜ P ) are indepen-dent of space and time. Hence we set ˜ ρ = ¯ ρ where¯ ρ = NL (53)is the average density. The action is minimised by tak-ing ˜ J = ( v L ¯ ρ/D ) ˜ P , note that this is O ( L ) when ex-pressed in these hydrodynamic variables, because theself-propulsion leads to ballistic motion. Then S = L (cid:20) γ D | f | + sω (¯ ρ, P ) D (cid:21) . (54)The next step is to use properties of ABPs to expressthe remaining quantities in terms of microscopic param-eters. Setting s = 0 in (54) and combining with (51),one can read off the probability of a trajectory where thepolarisation is fixed at P for all times between 0 and τ .This same probability can be computed applying large-deviation theory to the microscopic ABP model: the rel-evant LDP is given in (D8), and Appendix D 2 explainshow the rate function J can be computed as the Legen-dre transform of a Mathieu function, see also [35]. Hence S = L D [¯ ρ J ( P ) + sω (¯ ρ, P )] . (55)Recalling that ω is defined as the active work per unitvolume for a system with prescribed density and polarisa-tion, we have ω (¯ ρ, ) = ¯ ρ h w τ i . For P = 0, the consistentchoice is to define ω by considering an ensemble of trajec-tories that is biased by the polarisation. Details are givenin Appendix D 2, we define h·i h as an average analogousto (12), but with the bias acting on the polarisation, see(D10). Then ω (¯ ρ, P ) = h ¯ ρw τ i h ( P ) (56)where h ( P ) is defined by h ν τ i h ( P ) = P . For small P ,we show in Appendix D 4 that ω (¯ ρ, P ) = ¯ ρ h h w τ i + c ω | P | + O ( | P | ) i (57) c ω = ¯ ρτ L D r (cid:0) w τ , | ν τ | (cid:1) (58)where Cov( x, y ) = h xy i − h x ih y i is the covariance in thenatural (unbiased) ABP dynamics. Fig. 9(a) shows thatthis covariance is positive. That is, trajectories withlarger polarisation also tend to have larger active work.To analyse spontaneous symmetry breaking, we requirea Taylor expansion of S for small P . The behavior of J at small polarisation can be obtained exactly, the resultis (D18), which implies that J ( P ) = 12 D r | P | + O ( | P | ) . (59) . . . . . . . q . . . . . . . . . S s ( q ) s/D r =0 . − . − . − . − . s/D r =0 . − . − . − . − . r/σ C u , s ( ∆ r ) [ % ] − . s/D r . . . χ u , s (a)(b) FIG. 10. (a) Structure factor as function of wavevector k =2 π/λ , zoomed on highest computed wavelengths, at biasingparameter s . (b, main plot) Orientation correlation as afunction of the distance at biasing parameter s . (b, inset) Sum of the correlation function. Parameter values : N = 10 , l p /σ = 5, φ = 0 . n c = 10 , t max = 10 . Hence (55) becomes S = ¯ ρL D (cid:20) s h w τ i + 12 | P | ( D r + sc ω ) + O ( | P | ) (cid:21) . (60)Finally, minimising the action, we predict spontaneoussymmetry breaking for s < − s ∗ with s ∗ = D r /c ω . FromFig. 9(a), one sees that c ω depends weakly on D r , sothis is consistent with the observation of Sec. III, that s ∗ ∝ D r .Physically c ω > P , and so isthe rate function J associated with symmetry-breakingevents. Hence both contributions in (60) have the samescaling with P (and with L ), leading to a critical point s ∗ (independent of L ), for which the coefficients of thetwo relevant terms balance each other.Since this result was obtained by a top-down coarse-grained approach, one can expect that it should be quitegeneric, in that the same theory would be a natural de-scription of other active particles like RTPs. The keyingredient here is that the coefficient c ω in (57) shouldbe positive, which means by (58) that a larger polarisa-tion is correlated (in the unbiased steady state) with alarger active work. In ABPs, there is a clear mechanismfor this, that particles’ relative velocities are reduced iftheir align, because u i is a unit vector. In other activesystems, this would have to be checked on a case-by-casebasis.4 C. Fluctuations This mean-field analysis predicts that symmetry isspontaneously broken for s < − s ∗ with s ∗ > 0. Thismeans that the system remains isotropic in a finite rangearound s = 0. However, the fact that s enters (52) as sL indicates that the system can respond strongly to the biasalready for very small s , via hydrodynamic fluctuations.As explained in Appendix E, the (fast) polarisation fieldis not relevant on the hydrodynamic scale so it is suf-ficient for small bias to consider a scalar theory for thedensity. In this case density fluctuations and fluctuationsof the active work can be understood based on previousstudies [55, 71].The behaviour depends on the sign of ∂ ∂ρ ω (¯ ρ, ) whichwe abbreviate here by ¯ ω ′′ . For ABPs then ¯ ω ′′ < 0. Thefirst consequence of the hydrodynamic theory is that thevariance of w τ (under the natural dynamics) has a hy-drodynamic contribution. Appendix E shows that theasymptotic variance as L → ∞ behaves as¯ ρL lim τ →∞ τ Var( w τ ) = (¯ ω ′′ σ (¯ ρ )) π ¯ ρD c (¯ ρ ) [log L + O (1)] . (61)Fig. 9(b) shows numerical data that are consistent withthis prediction.The hydrodynamic theory of Appendix E, also predictsfor ω ′′ < s > s c with [42,55, 71] s c L = − π D c (¯ ρ ) ¯ ω ′′ σ (¯ ρ ) . (62)Recall that ω ′′ < s c > s , one expects the system tobe hyperuniform [55, 71]. Write ρ q for a Fourier com-ponent of the density field and S s ( q ) = h ρ q ρ − q i s for thestructure factor in the biased ensemble, where q = | q | .The structure factor is derived in (E16) which shows thatit behaves for small q as S s ( q ) = h ρ q ρ − q i s ≃ ( χ , s = 0 ,b s q, s < . (63)where χ = σ (¯ ρ ) /D c (¯ ρ ) and b s are constants [the behav-ior of b s and can be read from (E16)].Comparison of these theoretical predictions with nu-merical results requires some care because the numericalresults are limited to small systems, while the hydrody-namic theory is valid only if the system size is much largerthan all microscopic length scales.The existence of inhomogeneous states for s > s c wasalready discussed in [37]. For s ≤ 0, Fig. 10(a) showsthe behaviour of the structure factor in the isotropicphase. For s = 0, the system is a homogeneous fluid solim q → S ( q ) must be a positive constant. However, thenumerical results show S ( q ) increasing as q is reduced towards zero. The reason is that the system is not largeenough to show the hydrodynamic behavior – in a muchlarger system then smaller wavevectors would be acces-sible, and convergence of S ( q ) to its small- q limit wouldbe apparent. Similarly, for s < 0, one expects in verylarge systems to observe S ( q ) ∝ q at small q , but thisis not apparent from our numerical results in small sys-tems. Despite these restrictions, the numerical resultsshow that density fluctuations are suppressed in biasedensemble with s < 0, which is qualitatively consistentwith the theory, and with the physical reasoning that re-duced density fluctuations tend to suppress collisions andenhance the active work.Density fluctuations also have consequences for thesymmetry breaking transition at s = − s ∗ . First, thetransition takes place between a hyperuniform isotropicstate and a symmetry-broken state (which may also behyperuniform). In this case the simple relationship (48)should be adjusted, to account for the fact that the den-sity fluctuations near s ∗ are different from those of thebiased ensemble (D10). This effect presumably shifts thevalue of s ∗ but we expect qualitative predictions of mean-field theory to remain valid.An additional question is how density fluctuations cou-ple with those of the polarisation, close to the criticalpoint ( s = − s ∗ ). If one ignores the coupling between P and ρ in (52), one expects a transition in the universalityclass of an XY model in (2 + 1) = 3 dimensions (twospatial dimensions and one of time). This situation isfamiliar in quantum phase transitions for rotors [73].However, the coupling to a hydrodynamic density fieldmay change this universality class. For example, in theVicsek model, a locally-aligning interaction leads to clus-tering of the particles. As a result, symmetry breakingfor the orientations is coupled with collective motion andwith phase separation [14, 16]. By contrast, for the phasetransition considered here, clustering is suppressed (hy-peruniformity), see Fig. 10(a). Also, Fig. 10 shows thespatial orientation correlation function C u ,s (∆ r ) = R τ DP Ni,j =1 u ( θ i ( t )) · u ( θ j ( t )) δ ij (∆ r, t ) E s d t R τ DP Ni,j =1 δ ij (∆ r, t ) E s d tδ ij (∆ r, t ) = δ ( | r i ( t ) − r j ( t ) | − ∆ r ) (64)This average is performed over trajectories extractedfrom the cloning algorithm [49]. In these (small) systems,the main effect of the bias on fluctuations in the isotropicphase seems to be a constant (infinite-ranged) additivecontribution to C , see Fig. 10(b). This contribution is oforder 1 /N so that χ u ,s = Z C u ,s ( r ) 2 πr d r = 1 N * N X i =1 u ( θ i ( t )) ! + t,s (65)5is of order unity – and exactly equal to 1 in natural dy-namics ( s = 0). Such behaviour can be accounted for inmean-field theory – one might expect non-trivial spatialstructure to emerge in larger systems but analysis of sucheffects is beyond the scope of this work. VII. CONCLUSION We have analysed large deviations with enhanced ac-tive work in an ABP system. As found in [37], thisresults in spontaneous breaking of rotational symmetryand collective motion. We presented numerical evidenceand theoretical arguments that this transition occurs at s = − s ∗ with s ∗ ≃ D r , at least when D r is small (˜ l p ≫ w ∗ for theactive work, with w ∗ > h w τ i . This resolves an open ques-tion from [37], as to whether w ∗ = h w τ i .We have compared the behaviour of the CM statewith that of a controlled system where the ABP orien-tations interact via an infinite-ranged (mean-field) cou-pling. This controlled system captures the large devia-tions semi-quantitatively. Based on a hydrodynamic the-ory, we have also explained that we expect hyperuniformbehaviour for values of the active work between h w τ i and w ∗ . We discussed the extent to which the predictions ofthis theory should be generic in systems of self-propelledparticles.Despite this progress, several questions remain open,including the nature of the critical point where the sym-metry is broken. A related point is whether a controlledsystem with a more complicated (distance-dependent)coupling of orientations would capture the CM phasemore accurately. Another point of comparison is col-lective motion in aligning particles, such as the Vicsekmodel and its variants [14, 70, 74]. In that case, break-ing of rotational symmetry is often accompanied by phaseseparation into dense (polar) and more dilute (apolar) re-gions, which leads to a first-order transition to CM [16].For the systems considered here, the evidence [for ex-ample Fig. 10(a)] is that density fluctuations are sup-pressed in the CM phase, contrary to the enhanced den-sity fluctuations that one might expect in systems withpolar clusters. However, these numerical results for smallsystems are not sufficient to settle the nature of the phasetransition and its fluctuations. A detailed analysis ofsymmetry breaking within the Landau-Ginzburg theoryof Sec. VI would be an interesting direction for futurework. ACKNOWLEDGMENTS The authors acknowledge insightful discussions withTakahiro Nemoto, Julien Tailleur, and Thibaut Arnoulxde Pirey. This work was funded in part by the Euro-pean Research Council under the EU’s Horizon 2020 Pro-gramme, grant number 740269. ´EF acknowledges sup- port from an Oppenheimer Research Fellowship from theUniversity of Cambridge, and a Junior Research Fellow-ship from St Catharine’s College. MEC is funded by theRoyal Society.All the codes that were developed and used in thisproject are made freely available under the MIT license[75]. Appendix A: Large deviations for ABPs1. Eigenvalue problem Large deviations of w τ can be analysed through theeigenvalues of an operator called the backwards gener-ator. For compactness of notation we define (only forthis section) ˜ s = ( s/v ). Using results of [65, 76], theeigenvalue equation is ψ ( s ) F s = D X i ( ∇ i − ˜ s u i ) · ( ∇ i − ˜ s u i ) F s + X i ( v u i − D ∇ i U ) · ( ∇ i − ˜ s u i ) F s + D r X i ∂ F s ∂θ i (A1)where F s = F s ( r , . . . , r N , θ , . . . , θ N ) is the eigenvectorand ψ ( s ) the eigenvalue. This may be simplified as ψ ( s ) F s = X i [( v − sD ) u i − D ∇ i U ] · ∇ i F s + X i (cid:20) D ∇ i F s + D r ∂ F s ∂θ i (cid:21) + X i (cid:2) ˜ sD ( u i · ∇ i U ) + D ˜ s − ˜ s/v (cid:3) F s (A2)where the first two lines correspond to the generator ofan ABP system with a modified swim velocity v con s asdefined in (22), and the last line corresponds to a bias inwhich the only non-trivial term is v D ( u i · ∇ i U ), whichis the scalar product between the swim velocity and theinterparticle force, as it appears in the w f part of theactive work.The operator on the right hand side of (A2) is equal (upto an additive constant) to the operator that generatesthe biased ensemble on the right hand side of (21). Theadditive constant affects the eigenvalue but it does notaffect the biased ensemble. This is sufficient to establish(21).From the solution to the eigenproblem, one mayconstruct a corresponding optimally-controlled system,whose natural dynamics matches that of the biased en-semble. The corresponding optimal control potential is U opt s = − F s . (A3)6To construct the controlled system, consider the con-trolled dynamics (33). The backwards generator forthis model is L con , which acts on functions G = G ( r , . . . , r N , θ , . . . , θ N ) as L con ( G ) = X i (cid:20) D ∇ i G + D r ∂ G ∂θ i − D r ∂U con ∂θ i ∂ G ∂θ i (cid:21) + X i [ v con u i − D ∇ i ( U + U con )] · ∇ i G (A4)Denoting the right hand side of (A2) by L s ( F s ), thegenerator of the optimally-controlled dynamics may bethen be derived using the general formula L con ( G ) = F − s L s ( F s G ) − ψ ( s ) G , see for example [65]. Combiningthese ingredients, the optimally-controlled dynamics isgiven by (33,34). 2. KL divergence between controlled and naturaldynamics To derive a formula for the KL divergence in (39), werecall its definition D KL ( Q || P ) = Z Q ( x ) log Q ( x ) P ( x ) dx (A5)where Q, P are probability densities. Using standardtechniques in stochastic processes [77–79], the path prob-ability distribution for ABPs can be derived. We work inStratonovich calculus throughout. The path probabilityis P ( X ) ∝ exp[ − P i R τ S i ( X ( t )) dt ] where X = { r i , θ i } τ indicates a path and S i = | ˙ r i − v u i + D ∇ i U | D + ˙ θ i D r − D ∇ i U (A6)[We have neglected here a contribution to P ( X ) from theinitial condition, this does not cause a problem becausewe consider finally the long-time limit in (39).]Defining the corresponding quantity S con i for the con-trolled process one has D KL ( P con || P ) = N τ hS i − S con i i con (A7)where the average is computed in the steady state of thecontrolled process. One finds S i − S con i = v con − v D u i · ˙ r i + ∆ con i − | ( v con − v ) u i − D ∇ i U con | D + D ∇ U con − D r (cid:18) ∂U con ∂θ i (cid:19) + D r (cid:18) ∂ U con ∂θ i (cid:19) (A8)where ∆ con i has the property that P i R τ ∆ con i dτ =[ U con (0) − U con ( τ )] / U con appears in both the action and inthe equation of motion can be used to simplify the for-mulae for D KL , as for example in Eq. (H3) of [37]. Suchsimplifications are not essential for this work, so we omitthem.For the controlled dynamics of (37,38), it is useful todenote by ϕ the angle between ν and the x -axis. Thiscan be obtained by considering | ν | e i ϕ = 1 N X j e i θ j . (A9)Then the KL divergence of (39) can be obtained from(A7,A8), using ∇ i U con = 0 and v con = v , aslim τ →∞ N τ D KL ( P con g || P ) = (cid:28) g I ,τ − g D r I ,τ (cid:29) con − gN (A10)where I ,τ = 1 τ Z τ | ν ( t ) | d t, (A11) I ,τ = 1 N τ Z τ | ν ( t ) | X i sin( θ i ( t ) − ϕ ( t )) d t . (A12)The steady-state averages of these integrals are indepen-dent of τ ; hence the right hand side of (A10) is also in-dependent of τ and no limit is required there.The KL divergence of (A10) was evaluated by numer-ical simulation of the controlled dynamics, to obtain theupper bounds in Fig. 7. (Analytic results are availablefor this quantity in the large- N limit but numerical com-parison with the large-deviation rate function requiresresults for finite N .) Appendix B: Cloning algorithm with modifieddynamics1. Modified dynamics The cloning algorithm is applied to ABPs similarlyto [37]. The number of clones is denoted by n c . Wetypically repeat each computation N runs times, we takea simple average of the results and use the standard errorfor an error bar. Convergence with respect to n c is dis-cussed below. The computational cost of our calculationsis controlled primarily by n c .Given that we want to work with n c as small as pos-sible, the low probabilities of large deviation events canbe the origin of large systematic errors [48, 49, 80]. Toimprove the convergence of the algorithm, we use a mod-ification (or control) of the dynamics, informed by thetypical behaviour of the system for large fluctuations ofthe biasing observable. In particular we take (33) with v con = v con s from (22) and U con = U g from (37). Thechoice of the parameter g will be discussed below.7Following the same steps as (A8), we obtain the re-sult of [37] that the path probability distributions for thecontrolled and original dynamics are related as P [ X ] exp( − sN τ ( w τ )] ∝ P con [ X ] exp( − sN τ w mod τ )(B1)where X denotes a trajectory and the modified activework obeys sw mod τ = s (cid:18) − sDv + w f,τ (cid:19) − g (cid:18) N − I ,τ + gD r I ,τ (cid:19) (B2)where w f,τ is the force part of the active work from (8),and the integrals I ,τ and I ,τ are defined in (A11,A12).Physically, this means that the biased ensemble (12)for the ABP model (2) can be formulated alternativelyas a biased ensemble for the controlled ABP model, witha modified bias sw mod . Hence, the cloning algorithm isvalid as a method for sampling large deviations for anychoice of the parameter g (including g = 0).The role of the parameter g is to improve the numer-ical efficiency, and hence to obtain accurate results withsmaller numbers of clones. Several methods have beenproposed for determining suitable values of such param-eters [49, 80, 81]. Here we adopt the following method,similar to [37], where the value is chosen “on-the-fly”within the algorithm. We note that for the controlleddynamics with potential U g , we have for large N (and g < D r ) that h| ν | i con = 1 N − gD − r , (B3)The derivation is given in Appendix D 1. We obtainedgood results from the cloning algorithm by inverting thisequation for g and taking g = D r − N t − h R t d t ′ | ν ( t ′ ) | i clo ! (B4)where h . . . i clo designates an average over the cloneswithin the algorithm. Fig. 11 shows examples of suchtorque parameters as functions of the simulation time.Note, this always gives g < D r , the natural dynamicsof the controlled system is always in the paramagneticphase. However, as the biased system moves into the fer-romagnetic phase we find that g gets close to D r , andthe controlled system approaches the mean-field criticalpoint where fluctuations of ν are large. These large fluc-tuations are helpful for the algorithm, in that they gen-erate a wide range of trajectories, from which the cloningpart of the algorithm can select those with large valuesof w τ . 2. Convergence with respect to the number ofclones As stated in Ref. [49] the accuracy of the cloning algo-rithm is limited by the number n c of copies of our system − t − . − . − . − . . . g N = 50 D r = 0 . − . − . − . − . − . s/D r = 0 − . − . − . − . − . s/D r = 0 − t − . − . − . − . . . g N =1020304050 s/D r = − . D r = 0 . D r = 0 . (a) (b) FIG. 11. Torque parameter g as a function of time using therelation (B4) and setting g ( t = 0) = 0, for different (a) bias-ing parameters s and (b) number of particles N . Parametervalues : φ = 0 . t max = 10 , n c = 10 . s/D r = − . l p = 2 5˜ l p = 2 5 ∆ e rr h w τ i s [ ‰ ] s/D r = − . s/D r = − . n c [ × ] − s/D r = − . s/D r = − . ∆ e rr h | ν | τ i s [ ‰ ] n c [ × ]010 (a) (b) FIG. 12. Relative error of the value active work (left) andorder parameter (right) to their value at n c = 10 (leftmostvalue), ∆ err h•i s = (cid:2) h•i s − h•i s ( n c = 10 ) (cid:3) / h•i s ( n c = 10 ). Parameter values : N = 50, φ = 0 . t max = 10 , N runs = 10. which we simultaneously evolve. Most of the results pre-sented in the main text where computed for n c = 10 which is significantly lower than what was used for ex-ample in Ref. [37] (see Appendix D of that work).We show in Fig. 12 relative errors on the values of theactive work and the order parameter, for two values of thepersistence length and for N = 50 particles, when varyingthe number of clones from n c = 10 to n c = 5 · , withrespect to their value for the former number of clones.We have that the behaviour of the relative error on boththese quantities, which is at most of the order of 1%,indicates that the qualitative conclusions of the main textare robust. Appendix C: Biased trajectories of two RTPs on aring We consider the system of two RTPs defined in Sec. IV.The SCGF of (29) is the largest eigenvalue ψ RTP ( λ ) thatsolves ψ RTP ( λ ) P λ = ( L − λ ˙ w RTP f I ) P λ (C1)where L is the Fokker-Planck operator acting on a proba-bility distribution vector P λ ≡ ( P ++ λ , P −− λ , P + − λ , P − + λ ).The vector P λ depends on the distance r between the8particles, which is measured clockwise around the circu-lar system, starting at particle 1. Hence 0 ≤ r ≤ L . Since P is a vector, the operator L is a 4 × P α α = X α ′ ,α ′ = ± L α α ,α ′ α ′ P α ′ α ′ (C2)where the matrix elements of L can be read off from theFokker-Planck equation that corresponds to (26), whichis ˙ P α α = v ( α − α ) ∂∂r P α α + 2 ∂∂r (cid:18) P α α ∂∂r V (cid:19) + τ − (cid:0) P α α + P α α − P α α (cid:1) (C3)in which α i = − α i . In the hard-core limit, we expect theprobability density functions to take the form [64, 82] P α α λ ( r ) = ε α α λ ( r ) + γ α α , l λ δ ( r ) + γ α α , r λ δ ( L − r ) . (C4)Here ε α α λ is a smooth function of r that describes theprobability to find the particles with the given orienta-tions, at a separation r . Also γ α α , l λ and γ α α , r λ indicatethe probability that the particles are touching with par-ticle 1 to the left (l) or right (r) of particle 2, with theprescribed orientations.Particles never remain touching if their velocities pointaway from each other so γ + − , r λ = γ − + , l λ = 0. Note how-ever that particles may be touching but moving parallelto each other, γ ++ , r λ > w RTP f is non-zero only whenparticles are touching and have opposite orientation, so˙ w RTP f P + − λ = − v γ + − , l λ δ ( r )˙ w RTP f P − + λ = − v γ − + , r λ δ ( L − r ) (C5)with ˙ w RTP f P ++ λ = 0 = ˙ w RTP f P −− λ .The dominant eigenvector for (C1) obeys symmetryrelations stemming from particle interchangeability (C6)and parity symmetry (C7): P α α λ ( r ) = P α α λ ( L − r ) (C6) P α α λ ( r ) = P α α λ ( L − r ) . (C7)We take as normalisation condition Z L X α ,α = ± P α α λ ( r ) d r = 1 . (C8)The symmetry relations (C6, C7) imply for the prob-ability density function of aligned particles that ε ++ λ ( r ) = ε −− λ ( r ) = ε ααλ ( r ) (C9) ε ααλ ( r ) = ε ααλ ( L − r ) . (C10)Similarly for the “sticking” terms γ + − , r λ = γ − + , l λ = γ ααλ (C11) γ αα, l λ = γ αα, r λ = γ ααλ (C12) which are to be interpreted as definitions of γ ααλ , γ ααλ .These symmetry relations greatly simplify our calcula-tions.We highlight that the introduction of thermal diffusionin Eq. 26 would smooth out the δ -functions associated tosticking in Eq. C4, which will be replaced by exponentialsdecaying away from contact on a length scale set by thediffusivity [83]. 1. Unbiased steady state distribution We first solve the case without bias ( λ = 0) in steadystate ( ˙ P = 0). Evaluating Eq. C3 for r = 0 , L , gives0 = 2 v ∂∂r ε αα ( r ) + τ − [2 ε αα ( r ) − ε αα ( r )] . (C13)Using the symmetry properties (C6, C7) leads to ∂ r ε αα ( r ) = 0 and hence ε αα ( r ) = ε αα ( r ) = ε . In otherwords, P is independent of r and of the orientations, ex-cept when the particles are touching.Relations between γ αα and γ αα follow from integratingEq. C3 from r = 0 − to r = ǫ and then taking the hard-core limit: τ − ( γ αα − γ αα ) = 0 (C14) − v ε + τ − γ αα = 0 . (C15)Finally using the normalisation condition (C8)4 Lε + 4 γ αα + 2 γ αα = 1 . (C16)One can then solve to obtain P ++0 ( r ) = P −− ( r ) = a + laδ ( r ) + laδ ( L − r ) (C17) P + − ( r ) = a + 2 laδ ( r ) (C18) P − +0 ( r ) = a + 2 laδ ( L − r ) (C19)where a = [4( L + 2 l )] − , and the persistence length l = v τ p , as in the main text. This exactly recov-ers the continuous space and time results of Ref. [64](Eqs. 11, 12). Also (cid:10) ˙ w RTP f (cid:11) = − lv L + 2 l (C20)by (C5).Note that a steady state distribution that is non-uniform also for particles not in contact, as observed forRTPs on a discrete lattice [64], can be introduced byconsidering finite-time tumbles [84] or thermal diffusion[83]. 2. Biased steady state distribution We now solve the general biased case ( λ ∈ R ) using(C1). The method follows the unbiased case. For parti-9cles not in contact we have ψ RTP ( λ ) ε ααλ ( r ) = τ − ( ε ααλ ( r ) + ε ααλ ( L − r ) − ε ααλ ( r )) (C21) ψ RTP ( λ ) ε ααλ ( r ) = 2 αv ∂∂r ε ααλ ( r )+ τ − (2 ε ααλ ( r ) − ε ααλ ( r )) (C22)Hence d d r ( ε ααλ + ε ααλ ) = k λ ( ε ααλ + ε ααλ ) (C23)with k λ l = τ p ψ RTP ( λ )4 (cid:2) τ p ψ RTP ( λ ) (cid:3) . (C24)Note that k λ may be either real or imaginary. The solu-tions for ε (in both cases) can then be expressed as ε ααλ ( r ) = 12 + τ p ψ RTP ( λ ) A λ ( e − k λ r + e − k λ ( L − r ) ) (C25)2 ε ααλ ( r ) = (cid:18) − αk λ l τ p ψ RTP ( λ ) (cid:19) A λ e − k λ r + (cid:18) αk λ l τ p ψ RTP ( λ ) (cid:19) A λ e − k λ ( L − r ) (C26)where there is a single constant of integration A λ becausewe have enforced the symmetry (C10).We now derive four equations that can be used to ex-press ( γ ααλ , γ ααλ , A λ , λ ) as functions of ψ RTP ( λ ), whichenables a full solution of this problem.Integrating Eq. C1 from 0 − to ǫ , as in the unbiasedcase, and taking the ++ component of the vector P , oneobtains ψ RTP ( λ ) γ ααλ = τ − ( γ ααλ − γ ααλ ) . (C27)Similarly, taking the P − + component gives0 = − v ε − + λ (0 + ) + τ − γ ααλ , (C28)in which ε − + λ (0 + ) may be substituted using (C26) to ob-tain γ ααλ = τ p v A λ h (1+ e − k λ L )+ 2 k λ l τ p ψ RTP ( λ ) (1 − e − k λ L ) i . (C29)In addition, using (C25,C26) in the normalisation condi-tion (C8) leads to1 = 2 γ ααλ (cid:2) τ p ψ RTP ( λ ) (cid:3) + 2 A λ k λ (1 − e − k λ L )[4 + τ p ψ RTP ( λ )]2 + τ p ψ RTP ( λ ) (C30)where we have also used Eq. C27 to eliminate γ αα . Now, since the Fokker-Planck equation (C3) preservesthe normalisation of P , one may integrate (C1) over r and sum over α , α , then apply (C8) to obtain ψ RTP ( λ ) = 4 λv γ ααλ . Then use (C27) to obtain ψ RTP ( λ ) = 4 γ ααλ λv (cid:2) τ p ψ RTP ( λ ) (cid:3) . (C31)Eqs. (C27,C29,C30,C31) are the promised four equa-tions for ( γ ααλ , γ ααλ , A λ , λ ), in terms of ψ RTP ( λ ). Theproblem is now solved by computing the inverse of ψ RTP ( λ ), which amounts to treating ψ RTP as a parame-ter and solving for λ and the other variables. Note that k λ is fully determined by the value of ψ RTP , see (C24).To simplify the computation, we introduce dimension-less variables that treat the persistence length l as theunit of length: ˜ λ = λlv ˜ A ˜ λ = A λ l ˜ ψ RTP = τ p ψ RTP (C32)˜ k λ = k λ l ˜ L = Ll The ratio γ αα ˜ λ / ˜ A ˜ λ determines the relative probabilitiesof the particles being in contact or separated. One hasfrom (C29) that this ratio can be expressed in terms of˜ ψ RTP as γ αα ˜ λ ˜ A ˜ λ = 1 + e − ˜ k ˜ λ ˜ L k ˜ λ (1 − e − ˜ k ˜ λ ˜ L )2 + ˜ ψ RTP (C33)and, dividing Eq. C30 by Eq. C31 yields (after some re-arrangements):˜ λ = ˜ ψ RTP ( ˜ ψ RTP + 4)2( ˜ ψ RTP + 2) (cid:20) A ˜ λ γ αα ˜ λ (1 − e − ˜ k λ ˜ L )˜ k ˜ λ ( ˜ ψ RTP + 2) (cid:21) (C34)yielding ˜ ψ RTP (˜ λ ) ∼ λ when λ → ∞ . Combining(C33, C34) gives the promised inverse of ˜ ψ RTP (˜ λ ), as˜ λ = ˜ ψ RTP ( ˜ ψ RTP + 4)( ˜ ψ RTP + 2) " 12 + 1Ω ˜ L ( ˜ ψ RTP ) (C35)withΩ ˜ L ( ˜ ψ RTP ) = ˜ ψ RTP ψ RTP + 4)+ ˜ k ˜ λ ( ˜ ψ RTP + 2) 1 + e − ˜ k ˜ λ ˜ L − e − ˜ k ˜ λ ˜ L , (C36)where we used also (C24).Recall from (C24) that ˜ k is an imaginary number for − < ˜ ψ RTP < 0. In this case it is useful to rewrite˜ k ˜ λ e − ˜ k ˜ λ ˜ L − e − ˜ k ˜ λ ˜ L = | ˜ k ˜ λ | cot | ˜ k ˜ λ | ˜ L ! (C37)0With this result in hand, a careful analysis shows thatΩ ˜ L ( ψ ) has at least one zero for − < ψ < 0, at whichpoint ˜ λ diverges. This implies that ˜ ψ RTP (˜ λ ) has a hori-zontal tangent for ˜ λ → −∞ . The location of the (largest)zero sets the smallest possible value for ˜ ψ RTP (˜ λ ), whichis achieved as ˜ λ → −∞ . See Fig. 6(a). 3. Scaling regime It is instructive to consider two particles in a very largesystem, L → ∞ . The system has an associated scalinglimit whose behaviour is shown in Fig 13.Since the persistence length of the run-and-tumble mo-tion is much less than the system size, the particle motionon large scales can be characterised as (athermal) diffu-sion, and the particle explores the system on a time scale O ( L ). Since ψ RTP is an inverse time scale, it is expectedto be O ( L − ). Moreover, it follows from Eq. (C20) thattypical values of w RTP f are of order L − in this regime.Physically, this small value arises because of the smallfraction of time that the particles spend in contact, when L is large. Hence ψ RTP ( λ ) = 2 lv L + 2 l λ + O ( λ ) . (C38)One then expects a scaling form as L → ∞ : ψ RTP ( λ ) ≃ L − ϕ ( λL ) , (C39)which will be verified below. The corresponding form ofthe rate function for w RTP f is obtained from (14) as I ( w f ) ≃ L − I ( w f L ) . (C40)The natural dimensionless quantities in this regime areΛ = λLv τ p L l ψ RTP ( λ ) (C42)so that ϕ ( λL ) = 2 l Ψ(Λ) /τ p . The quantities Λ , Ψ(Λ)have same sign.Since ψ RTP is O ( L − ) then Eq. C35 shows that k = O ( L − ). At the lowest order in L − we then infer fromEq. C35 that for Λ > √ Ψ tanh (cid:16) √ Ψ (cid:17) , (C43)while for Λ < − p | Ψ | tan (cid:16)p | Ψ | (cid:17) (C44)yielding lim Λ →−∞ Ψ(Λ) = − π / L − , it follows fromEqs. C25, C26 that ε αα Λ ( r ) = ε αα Λ ( r ) = ε Λ ( r )= 12 A Λ (cid:16) e − k Λ r + e − k Λ ( L − r ) (cid:17) (C45) . . . . . . r/L . . . . . L ε Λ ( r ) Λ = − 3Λ = − − 1Λ = 0 Λ = 1Λ = 2 FIG. 13. Regular part of the probability density function ε Λ ( r ) from (C48,C49) scaled by the ring length L . and from Eqs. C31, C33 that γ Λ = l L Ψ(Λ)Λ (C46) A Λ = 2 γ Λ l 11 + e − k Λ L . (C47)Using Eq. C24 leads for Λ > ε Λ ( r ) = 14 L Ψ(Λ)Λ cosh (cid:16)p Ψ(Λ) (cid:0) − rL (cid:1)(cid:17) cosh (cid:16)p Ψ(Λ) (cid:17) (C48)and for Λ < ε Λ ( r ) = 14 L Ψ(Λ)Λ cos (cid:16)p | Ψ(Λ) | (cid:0) − rL (cid:1)(cid:17) cos (cid:16)p | Ψ(Λ) | (cid:17) (C49)which we plot in Fig. 13.The physical picture emerging from Fig. 13 is as fol-lows. In a large system, a very weak bias λ = O ( L − )is sufficient to change qualitatively the separation of theparticles. The resulting probability distributions are in-dependent of particle orientation but depend on the par-ticle separation through the scaling variable r/L . ForΛ > < L − , this holds throughout the scal-ing regime λ = O ( L − ). 4. Distribution over the infinite-time interval The probability density vector P λ which satisfies (C1,C8) indicates the fraction of trajectories for which theparticles have final orientations α i and separation r in the λ -ensemble [49, 65]. We are also interested in the fractionof time spent with given orientations α i and separation r in the λ -ensemble, which we will denote ˆ P λ . We ex-pect these probability density functions to take the same1form and respect the same symmetries as their final timecounterparts (C4, C6, C7), so thatˆ P ααλ ( r ) = ˆ ε ααλ ( r ) + ˆ γ ααλ δ ( r ) + ˆ γ ααλ δ ( L − r ) (C50)ˆ P + − λ ( r ) = ˆ ε + − ( r ) + ˆ γ αα δ ( r ) (C51)ˆ P − + λ ( r ) = ˆ ε − + ( r ) + ˆ γ αα δ ( L − r ) (C52)in addition to being normalised accordingly (C8). In or-der to compute them it is necessary to solve the eigen-problem adjoint to (C1) ψ RTP ( λ ) Q λ = ( L † − λ ˙ w RTP f I ) Q λ (C53)where the same ψ RTP ( λ ) is the largest eigenvalue. It thenfollows that ˆ P α α λ ( r ) = P α α λ ( r ) Q α α λ ( r ) (C54)according to Ref. [65].The matrix elements of L † can be read off from thebackward Fokker-Planck equation that corresponds to(C3)˙ Q α α = − v ( α − α ) ∂∂r Q α α − ∂∂r Q α α ∂∂r V + τ − (cid:0) Q α α + Q α α − Q α α (cid:1) (C55)and indicate that the eigenvector of this operator corre-sponding to the eigenvalue 0, which is also Q in (C53),is constant. Since bias introduces in (C53) terms whichare at least as regular as those in (C55), we expect Q α α λ to remain smooth and continuous for λ = 0. For parti-cles not in contact, Q α α λ satisfies the same equations as P α α λ (C21, C22), with the replacement v → − v . Wecan then finally concludeˆ ε ααλ ( r ) = ˆ A λ A λ ε ααλ ( r ) (C56)ˆ ε ααλ ( r ) = ˆ A λ A λ ε ααλ ( r ) ε ααλ ( r ) (C57)ˆ γ ααλ = ˆ A λ A λ γ ααλ ε ααλ (0 + ) (C58)ˆ γ ααλ = ˆ A λ A λ γ ααλ ε + − λ ( L − ) (C59) where ˆ A λ is a normalisation constant for ˆ P λ . 5. Polarisation From (31,32) we identify ν RTPave ( λ ) = Z L (cid:16) ˆ P ++ λ ( r ) + ˆ P −− λ ( r ) (cid:17) d r (C60)and (replacing ˆ P with P ) ν RTPend ( λ ) = Z L (cid:0) P ++ λ ( r ) + P −− λ ( r ) (cid:1) d r . (C61)An exact expression of the polarisation ν RTPend (C61) canbe computed by noting that ν RTPend ( λ ) = 4 γ ααλ + 2 R L ε ααλ ( r ) d r γ ααλ + 2 γ ααλ + 2 R L [ ε ααλ ( r ) + ε ααλ ( r )] d r (C62)where the denominator is 1 by (C8) and we used that R L [ ε ααλ ( r ) − ε ααλ ( r )] dr = 0. Moreover, Eq. C27 yields γ ααλ = ( ˜ ψ RTP + 2) γ ααλ and Eqs. C25, C26 yield2 Z L ε ααλ ( r ) d r = ( ˜ ψ RTP + 2) Z L ε ααλ ( r ) d r (C63)therefore ν RTPend = 2˜ ψ RTP + 4 (C64)whose dependence on ˜ λ can then be obtained parametri-cally via (C35).An exact expression of the polarisation ν RTPave (C60) isalso available by writing ν RTPave ( λ ) = 4 γ ααλ ε ααλ (0 + ) + 2 R L ε ααλ ( r ) d r γ ααλ ε ααλ (0 + ) + 2 γ ααλ ( ˜ ψ RTP + 2) ε + − λ ( L − ) + 2 R L ε ααλ ( r ) d r + 2 R L ε ααλ ( r ) ε ααλ ( r ) d r (C65)which avoids the need to determine ˆ A λ , and where we have used (C56, C57, C58, C59). This allows ν RTPave to bedetermined in full via (C25, C26, C33).It is easily verified that ν RTPend = ν RTPave = 1 / λ = 0, corresponding to a state where aligned and anti-aligned2states are equiprobable. Appendix D: Statistics of orientational orderparameter(s) This Appendix analyses (controlled) ABPs in situa-tions where their orientations evolve independently oftheir positions. In these cases the statistics of the or-der parameter ν can be computed. 1. Mean-field analysis We consider the dynamics of the particle orientationsalone, for the controlled system (37). We have U con g = − gN − P ij cos( θ i − θ j ) so the controlled equation of mo-tion for the orientation is˙ θ i = gN − X j sin( θ i − θ j ) + p D r ξ i (D1)Writing sin( θ i − θ j ) = sin θ i cos θ j − cos θ i sin θ j , and using(16) with ν = | ν | (cos ϕ, sin ϕ ) yields (38) of the maintext.For N ≫ 1, mean-field theory is valid, because theindividual θ i relax much faster than the global ϕ , andfluctuations of ν are also negligible. Without loss of gen-erality, we take ϕ = 0, hence it is consistent to set | ν | = h cos θ i i con (D2)in (38). Treating this quantity as a fixed number, thesteady state of the system obeys a Boltzmann distribu-tion where each θ i is independent with distribution p con | ν | ( θ i ) ∝ exp (cid:18) g | ν | D r cos θ i (cid:19) . (D3)where the constant of proportionality is fixed by normal-isation. Combining (D2,D3) leads to a self-consistencyrelationship | ν | = Z p con | ν | ( θ ) cos θ d θ (D4)The integral can be expressed in terms of a Bessel func-tion. However, the relevant question is for which valuesof g non-trivial solutions exist (excluding | ν | = 0). Forthat purpose one may expand for small values of the pa-rameter 2 g | ν | /D r which yields h cos θ i i con = ( g/D r ) | ν | + O ( | ν | ). The correction term is negative and hence thenon-trivial (ferromagnetic) solution appears for g > D r . (D5)To analyse the paramagnetic phase we have by thecentral theorem of Sec. III B that for g = 0 then p ( ν ) ∝ e − N | ν | . The Boltzmann distribution for the steady state of the controlled system is obtained by multiplying bye − U con g , yielding p g ( ν ) ∝ exp (cid:2) − N | ν | (cid:0) − gD − r (cid:1)(cid:3) (D6)from which we obtain (B3). Consistent with the previousargument, that fluctuation diverges at the critical point g = D r . For larger g , estimation of p by the central limittheorem is too simplistic and a more detailed analysis isrequired, for example as in (D2,D3). 2. Large deviations of the time-averaged orderparameter In this section we consider large deviations of the time-averaged (vectorial) order parameter ν τ = 1 τ Z τ ν ( t ) dt . (D7)Recall that ν ( t ) is defined in (16) as a simple averageof individual orientations which evolve independently by(2). Hence the statistics of ν τ can be analysed exactly.As τ → ∞ there is a LDP p ( ν τ ) ∼ exp[ − τ N J ( | ν τ | )] (D8)where J is the rate function, which only depends on themodulus of ν τ , by symmetry. Note that the function J isdistinct from J in (40), which describes large deviationsof the time-integrated modulus of ν . See however (D22),below. There is an associated SCGF ψ OP ( h ) = lim τ →∞ N τ log h exp( − τ N h · ν τ ) i (D9)where h = | h | ; the right hand side only depends on themodulus of h , by symmetry. There is a correspondingbiased ensemble of trajectories, in which a generic ob-servable A has average value hAi h = hA exp( − τ N h · ν τ ) ih exp( − τ N h · ν τ ) i . (D10)Since the rotors are independent under the ABP dy-namics, the expectation value in (D9) reduces to a prod-uct of expectation values for single rotors. Hence ψ OP solves the eigenvalue problem ψ OP ( h ) F h ( θ ) = D r F ′′ h ( θ ) − h F h ( θ ) cos θ (D11)As noted in [35], this problem is related to Mathieu’sequation. Let ˜ θ = θ/ h -dependent quantities a = − ψ OP ( h ) /D r and q = 2 h/D r . Defining also ˜ F (˜ θ ) = F h (2˜ θ ) we have˜ F ′′ (˜ θ ) + ( a − q cos 2˜ θ ) ˜ F (˜ θ ) = 0 (D12)We recognise Eq. D12 as Mathieu’s differential equation[85]. For any real number q there is a countable infinity3of possible values of a and associated solutions ˜ F . We areinterested in functions F h that are even and 2 π -periodic.Hence ˜ F must be π -periodic in ˜ θ . We therefore introduce M (0) (˜ θ, q ) which is the 0 th even and π -periodic Mathieufunction and a (0) M ( q ) its characteristic value [35], such that ψ OP ( h ) = − D r a (0) M (2 h/D r ) , (D13) F h ( θ ) = M (0) (cid:18) θ , hD r (cid:19) . (D14)Hence by Legendre transform the rate function in (D8)is J ( ν ) = sup h [ − hν − ψ OP ( h )] . (D15)This result is exact for all N and all ν .To obtain additional physical insight, we obtain thequadratic behaviour of the rate function at small ν .This requires that we solve (D12) for small q , whichis a computation in perturbation theory [86]. Since˜ F is even and π -periodic in ˜ θ it can be expanded as˜ F (˜ θ ) = 1 + P ∞ n =1 β n cos 2 n ˜ θ . The coefficients β n andthe eigenvalue a can then be expanded in powers of q .To leading order, a = − q / O ( q ) , β = − q/ O ( q ) (D16)and β n = O ( q n ) for n ≥ 2. Hence by (D13,D14) ψ OP ( h ) = 12 D r h + O ( h ) , (D17)and so for small ν , J ( ν ) = 12 D r ν + O ( ν ) . (D18)The corresponding eigenfunction is F h ( θ ) = 1 − hD r cos θ + O ( h ) , (D19)Hence the optimal control potential (for this single ori-entation vector) is U OPopt = − F h , so U OPopt ( θ ) = 2 hD r cos θ + O ( h ) . (D20)Recall, we are considering here large deviations wherethe order parameter is aligned parallel (or anti-parallel)the x -axis. For h > 0, the the control potential acts toalign the orientation vectors anti-parallel to this axis, asexpected from (D10). 3. Time-averaged modulus of the order parameter The discussion of the LDP (D8) of the previous sectionis useful as a way to characterise the LDP (40) of themain text. In this case the relevant SCGF is ψ ( λ ) = lim τ →∞ N τ log h exp( − τ N λν τ ) i (D21) . . . . . ν . . . . . . . . J ( ¯ ν ) N = 10 N = 20 N = 30 N = 40 N = 50 J (¯ ν ) D r = 0 . D r = 0 . . . . . . √ N ¯ ν . . . . . . N J ( ¯ ν ) N = 10 N = 20 N = 30 N = 40 N = 50 D r = 0 . D r = 0 . (a) (b) FIG. 14. (a, b) Rate function of the time-averaged modulusof the order parameter J (¯ ν ) computed from cloning simula-tions of rotors ( Parameter values : n c = 10 , t max = 10 ) (a) and compared to the rate function of the time-averaged orderparameter J (¯ ν ) (See Eq. D15). and J ( ν ) = sup λ [ − λν − ψ ( λ )]. In contrast to the pre-vious case this problem cannot (to our knowledge) besolved exactly for finite N . However as N → ∞ theproblem is of mean-field type. In this case, the intuitiveresult is that the large-deviation mechanism for ν τ shouldbe the same as that of ν τ , so thatlim N →∞ J ( ν ) = J ( ν ) (D22)as illustrated by Fig. 14(a).We note from this figure that J has its zero at theorigin, h ν τ i = 0. On the other hand, (18) shows that theunique zero of J is at h ν τ i = p π/ (4 N ) and in fact J increases rapidly for smaller values of ν τ . Hence, for anygiven N , there is a region between 0 and h ν τ i where J deviates strongly from J . However, this region vanishesas N → ∞ so (D22) holds for all ν > ψ ( k ) = lim τ →∞ τ log D e − k R τ √ N | ν ( t ) | d t E (D23)is a well-defined smooth function for N ≫ 1. Therefore,with N ψ ( λ ) = ˜ ψ ( √ N λ ), we obtain N J ( ν ) ≃ B (cid:16) √ N ν − p π/ (cid:17) (D24)in the regime where √ N ν = O (1), and where B is aconstant independent of N , as illustrated by Fig. 14(b).The next step is to outline a derivation of (D22),which also yields the optimally-controlled dynamics forthis problem. The SCGF can be obtained by solving aneigenproblem N ψ ( λ ) F pol λ = D r X i ∂ ∂θ i F pol λ − λN | ν |F pol λ (D25)where F pol λ = F pol λ ( θ , . . . , θ N ) depends on all angles.Since this is a mean-field problem the solution for large N can be approximated as F pol λ = f ( ν ) Y i ζ λ ( θ i | ν ) (D26)4where the orientation vectors interact only through theiraverage, which is ν . Note also | ν | = N − P i cos( θ i − ϕ ) where ϕ is the angle between ν and the x -axis. Weassume that ζ is normalised as R ζ ( θ | ν ) d θ = 1. Thefunction f is assumed to depend only on | ν | , and forlarge N it is sharply-peaked in this variable, at | ν | = ν ∗ .In order that F pol is also sharply peaked at ν ∗ , we requirea self-consistency condition Z cos( θ − ϕ ) ζ ( θ | ν ) d θ = ν ∗ (D27)The eigenproblem (D25) is Hermitian (self-adjoint)so there is a variational (Rayleigh-Ritz) formula for itslargest eigenvalue. Using F pol as ansatz yields ψ ( λ ) & − λν ∗ − πN X i Z Ψ ζ ( θ i , ν ∗ )d ϕ (D28)where we used (D27) as well as ν ∗ = ν ∗ (cos ϕ, sin ϕ ) andΨ ζ ( θ, ν ) = − Z D r ζ ( θ | ν ) ∂ ∂θ ζ ( θ | ν )d θ i . (D29)We have neglected terms in (D28) which arise from theaction of the derivatives on ν , these are negligible as N → ∞ . In fact, the mean-field structure of the problemmeans that the bound (D28) will become an equality as N → ∞ , if the optimal choice is made for ζ .It is convenient to work in terms of the rate function,using J ( ν ) = sup λ [ λν − ψ ( λ )] to see that J ( ν ) ≃ inf ζ Ψ ζ ( θ, ν ) where the maximisation is again subject to(D27). Implementing this constraint and the normalisa-tion constraint on ζ by Lagrange multipliers µ , µ , wefind D r ∂ ζ∂θ + [ µ − µ cos( θ − ϕ )] ζ = 0 (D30)Hence we have recovered Mathieu’s equation. Proceedingsimilar to Appendix D 2 and using that extremisation ofthe Lagrange multipler µ yields a maximimum in thiscase, one obtains J ( ν ) ≃ sup µ [ − µ ν − ψ OP ( µ )] . (D31)The notation ≃ indicates that this relation becomes exactas N → ∞ . Eq. (D22) follows on comparing with (D15).It follows that the optimal control potential for the LDPof (40,D21) is U polcon ( θ , . . . , θ N ) = − X i log F λ ( θ i − ϕ ) (D32)where the notation F λ indicates the function defined in(D14), evaluated at h = λ . This is indeed a mean-field-type interaction among orientations. By the same argu-ment as (D20), it reduces for small λ to U polcon ( θ , . . . , θ N ) ≃ − λN D r X ij cos( θ i − θ j ) . (D33)This is nothing but U con g from (37) with g = 2 λ/D r . Theresult is that for small values of ν , Eq. (37) is an optimalcontrol potential for large deviations of the orientation. 4. Expansion of ω We assume that ω (¯ ρ, P ) can be inferred from the en-semble of trajectories biased with respect to the polari-sation P (Eq. 56), ω (¯ ρ, P ) = (cid:10) ¯ ρw τ e − h ( P ) · τN ¯ ν τ (cid:11)(cid:10) e − h ( P ) · τN ¯ ν τ (cid:11) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h ( P ) , h ¯ ν τ i h ( P ) = P (D34)where we have used the averaged polarisation fromEq. D7.We have in the limit h ( P ) → , D w τ e − h ( P ) · τN ¯ ν τ E = h w τ i + 14 τ N | h ( P ) | (cid:10) w τ | ¯ ν τ | (cid:11) , (D35) D e − h ( P ) · τN ¯ ν τ E = 1 + 14 τ N | h ( P ) | (cid:10) | ¯ ν τ | (cid:11) (D36)up to O ( h ( P ) ) terms and where we have discarded lin-ear terms in h ( P ) by symmetry, therefore¯ ρ − ω (¯ ρ, P ) = h w τ i + 14 τ N | h ( P ) | Cov( w τ , | ¯ ν τ | )(D37)linking ω (¯ ρ, P ) and the covariance of the active work andthe squared average polarisation.We note that h ¯ ν τ i h ( P ) = − τ N h ( P ) Var( ¯ ν τ ) (D38)using h ¯ ν τ i = in the limit τ → ∞ such that Var( ¯ ν τ ) = (cid:10) | ¯ ν τ | (cid:11) , and (cid:10) | ¯ ν τ | (cid:11) = 2 D r τ N (cid:18) − D r τ (cid:0) − e − D r τ (cid:1)(cid:19) = 2 D r τ N , τ → ∞ (D39)using h u ( θ i ( t )) · u ( θ j ( t ′ )) i = δ ij e − D r | t − t ′ | from Eq. 2,therefore with h ¯ ν τ i h ( P ) = P we may write h ( P ) = − τ N Var( ¯ ν τ ) P = τ →∞ − D r P , (D40)linking the biasing parameter and the polarisation.We then have¯ ρ − ω (¯ ρ, P ) − h w τ i = | P | ν τ ) Cov( w τ , | ¯ ν τ | )= 14 | P | τ N D r Cov( w τ , | ¯ ν τ | ) , τ → ∞ , (D41)at leading order in | P | .5 Appendix E: Fluctuations of the active work in thehydrodynamic theory This Appendix describes density fluctuations of ABPsat hydrodynamic level, including large deviations. Wefollow Ref. [55], which draws on earlier results includ-ing [71, 72]. At hydrodynamic level, we are restricted tosmall biasing fields s = O (1 /L ). In this case, the (fast)polarisation field is unaffected by the bias and can besafely integrated out. At the level of (46), this leads torenormalisation of the diffusion constant D c but we donot distinguish the bare and renormalised values of D c ,for simplicity. Since the polarisation has been integratedaway, this analysis of density fluctuations is restricted tostates with h P i = 0, but this is sufficient to cover homo-geneous phases for s > s < 1. Quadratic theory As in Refs. [55] (Section 5.1 and Appendix B), we con-sider a perturbation around the homogeneous profile ρ ( r , t ) = ¯ ρ + δρ ( r , t ) (E1)with δρ ≪ ¯ ρ and R δρ ( r , t )d r = 0. Since we consider anisotropic system, P = , we define¯ ω = ω (¯ ρ, P = ) (E2)¯ ω ′′ = ∂ ∂ρ ω (¯ ρ, P = ) (E3)so that we can Taylor expand the active work over ρ , N τ w τ = L τ ¯ ω + 12 ¯ ω ′′ Z [0 ,L ] Z τ ( δρ ) d r d t + O ( δρ ) . (E4)At the consistent level of expansion, the stochastic equa-tion for the density (46, 47) is ∂∂t δρ = D c (¯ ρ ) ∇ δρ − p σ (¯ ρ ) ∇ · η . (E5)We introduce the Fourier modes of the density,˜ ρ q = 1 L Z [0 ,L ] δρ ( r ) e − i q · r d r (E6)so that δρ = X q =(0 , ˜ ρ q e i q · r . (E7)Hence N τ w τ = L τ ¯ ω + L ¯ ω ′′ X q x ≥ ,q y q =(0 , Z τ ˜ ρ q ˜ ρ − q d t (E8) where the sum runs over non-zero modes q =2 π ( n x , n y ) /L , with q x ≥ 0, and where we have used Z [0 ,L ] δρ d r = L X q =(0 , ˜ ρ q ˜ ρ − q = 2 L X q x ≥ ,q y q =(0 , ˜ ρ q ˜ ρ − q (E9)according to Parseval’s theorem. Since the theory is de-fined on the mesoscopic scale, sums over q are restrictedto | q | < Λ where Λ is an upper cutoff of order unity [tobe precise, it is of order | Ω r | − /d , for consistency with(44)].From Eqs. E5 and E6, we derive the stochastic equa-tion satisfied by the non-zero Fourier modes ∂∂t ˜ ρ q = − D (¯ ρ ) | q | ˜ ρ q + p σ (¯ ρ ) | q | ˜ η q (E10)where the longitudinal part of the noise term, ˜ η q , is acomplex Gaussian white noise with zero mean and vari-ance h ˜ η q ( t )˜ η ∗ q ( t ′ ) i = 1 L | q | × Z Z [0 ,L ] h [( − i q ) · η ( t )][(i q ) · η ( t ′ )] i d r d r ′ = δ ( t − t ′ ) , (E11)and also hℜ (˜ η q ( t )) ℑ (˜ η q ( t ′ )) i = 0 . Noises with differentwavevectors q are independent. The key point is that(E10) is diagonal in q so every wavevector can be anal-ysed separately. 2. Biased ensemble of trajectories The next step is to consider a biased ensemble de-fined by the (linear) equations of motion (E10) and thereweighting factor e − sτNw τ , where the exponential factor(E8) is quadratic (and diagonal) in the density fluctua-tions.Following Ref. [55] (Appendix B), we consider the com-plex Ornstein-Uhlenbeck process z = − ζz + p γη (E12)where η is a complex Gaussian white noise withthe same statistics as ˜ η q . For a biased ensemblewith exponential biasing factor of e − sK τ where K τ = α R τ | z ( t ) | dt , the scaled cumulant generating functionfor K τ can be computed. Identifying ( z, ζ, γ, α ) with( ρ q , D c (¯ ρ ) q , σ (¯ ρ ) q , ω ′′ ), the resulting SCGF for w τ isobtained by summing over the modes, to find6 ψ ( s ) = lim τ →∞ N τ log (cid:10) e − sNτw τ (cid:11) = − s h w τ i − N X q x ≥ ,q y q =(0 , (cid:18)q D c (¯ ρ ) | q | + 2 s ¯ ω ′′ σ (¯ ρ ) | q | − D c (¯ ρ ) | q | − s ¯ ω ′′ σ (¯ ρ ) D c (¯ ρ ) (cid:19) (E13)consistent with [71, 72]. We emphasise that this result isvalid only on the hydrodynamic scale, which means verysmall bias s = O (1 /L ).Several results are available from this formula. We firstcompute − w ′ ( s ) = ψ ′′ ( s ) = [ ω ′′ σ (¯ ρ )] D (¯ ρ ) N X q x ≥ ,q y q =(0 , q (E14)which is related to the variance of the active work fromEq. 20. The sum in this last expression can be approxi-mated as X q x ≥ ,q y q =(0 , q ≃ L (2 π ) Z Λ2 π/L πq d qq = L π [log L + O (1)] (E15)where Λ is the upper cutoff on q . Hence (61) follows,by combining these results with (20,53). The origin of this diverging variance is the presence of a slow (hydro-dynamic) time scale, diverging as L .The second result that is available from (E13) is thatthe argument of the square root will become negative if sω ′′ is sufficiently negative, indicating that the systembecomes inhomogeneous. The instability is in the low-est mode, which has | q | = 2 π/L . Noting that ω ′′ < ψ in (E13) is the scaledcumulant generating function for squared density fluctu-ations, the structure factor of the biased ensemble canalso be obtained by taking a derivative, leading to h| ρ q | i s = σ (¯ ρ ) | q | p D c (¯ ρ ) | q | + 2 s ¯ ω ′′ σ (¯ ρ ) | q | (E16)similar to [71]. Observe that the limiting behaviour ofthis function at q → sω ′′ is zero or positive. In the latter case then h| ρ q | i s → q → 0, corresponding to hyperuniformity. Hence(63) follows. If sω ′′ < q → q , which signals phase separation,as noted above. [1] M. C. Marchetti, J. F. Joanny, S. Ramaswamy,T. B. Liverpool, J. Prost, Madan Rao, andR. Aditi Simha, “Hydrodynamics of soft active matter,”Rev. Mod. Phys. , 1143–1189 (2013).[2] Clemens Bechinger, Roberto Di Leonardo, HartmutL¨owen, Charles Reichhardt, Giorgio Volpe, and Gio-vanni Volpe, “Active particles in complex and crowdedenvironments,” Rev. Mod. Phys. , 045006 (2016).[3] ´Etienne Fodor and M. C.Marchetti, “The statisticalphysics of active matter: From self-catalytic colloids toliving cells,” Physica A , 106 – 120 (2018).[4] Xiao-Lun Wu and Albert Libchaber, “Particle dif-fusion in a quasi-two-dimensional bacterial bath,”Phys. Rev. Lett. , 3017–3020 (2000).[5] J Elgeti, R G Winkler, and G Gompper, “Physics ofmicroswimmerssingle particle motion and collective be-havior: a review,” Rep. Prog. Phys. , 056601 (2015).[6] Andrea Cavagna, Alessio Cimarelli, Irene Gi-ardina, Giorgio Parisi, Raffaele Santagati,Fabio Stefanini, and Massimiliano Viale,“Scale-free correlations in starling flocks,”Proc. Natl. Acad. Sci. USA , 11865–11870 (2010).[7] Andrea Cavagna and Irene Giardina, “Bird flocks as condensed matter,”Ann. Rev. Condens. Matter Phys. , 183–207 (2014).[8] Arianna Bottinelli, David T. J. Sumpter, and Jesse L.Silverberg, “Emergent structural mechanisms for high-density collective motion inspired by human crowds,”Phys. Rev. Lett. , 228301 (2016).[9] Nicolas Bain and Denis Bartolo, “Dynamic re-sponse and hydrodynamics of polarized crowds,”Science , 46–49 (2019).[10] Julien Deseigne, Olivier Dauchot, and HuguesChat´e, “Collective motion of vibrated polar disks,”Phys. Rev. Lett. , 098001 (2010).[11] Nitin Kumar, Harsh Soni, Sriram Ramaswamy, andA. K. Sood, “Flocking at a distance in active granularmatter,” Nat. Commun. , 4688 (2014).[12] Jeremie Palacci, Stefano Sacanna, Asher Preska Stein-berg, David J. Pine, and Paul M. Chaikin,“Living crystals of light-activated colloidal surfers,”Science , 936–940 (2013).[13] Ivo Buttinoni, Julian Bialk´e, Felix K¨ummel, Hart-mut L¨owen, Clemens Bechinger, and ThomasSpeck, “Dynamical clustering and phase separationin suspensions of self-propelled colloidal particles,” Phys. Rev. Lett. , 238301 (2013).[14] Tam´as Vicsek, Andr´as Czir´ok, Eshel Ben-Jacob,Inon Cohen, and Ofer Shochet, “Novel type ofphase transition in a system of self-driven particles,”Phys. Rev. Lett. , 1226–1229 (1995).[15] John Toner and Yuhai Tu, “Long-range order in a two-dimensional dynamical XY model: How birds fly to-gether,” Phys. Rev. Lett. , 4326–4329 (1995).[16] H. Chat´e, “Dry aligning dilute active matter,”Annu. Rev. Condens. Matter Phys. , 189–212 (2020).[17] Yaouen Fily and M. Cristina Marchetti, “Athermal phaseseparation of self-propelled particles with no alignment,”Phys. Rev. Lett. , 235702 (2012).[18] Gabriel S. Redner, Michael F. Hagan, andAparna Baskaran, “Structure and dynamicsof a phase-separating active colloidal fluid,”Phys. Rev. Lett. , 055701 (2013).[19] Michael E Cates and Julien Tailleur,“Motility-induced phase separation,”Annu. Rev. Condens. Matter Phys. , 219–244 (2015).[20] Xingbo Yang, M. Lisa Manning, and M. CristinaMarchetti, “Aggregation and segregation of confined ac-tive particles,” Soft Matter , 6477–6484 (2014).[21] S. C. Takatori, W. Yan, and J. F. Brady,“Swim pressure: Stress generation in active matter,”Phys. Rev. Lett. , 028103 (2014).[22] Alexandre P. Solon, Joakim Stenhammar, RaphaelWittkowski, Mehran Kardar, Yariv Kafri, Michael E.Cates, and Julien Tailleur, “Pressure and phaseequilibria in interacting active brownian spheres,”Phys. Rev. Lett. , 198301 (2015).[23] ´Etienne Fodor, Cesare Nardini, Michael E. Cates,Julien Tailleur, Paolo Visco, and Fr´ed´eric van Wi-jland, “How far from equilibrium is active matter?”Phys. Rev. Lett. , 038103 (2016).[24] Dibyendu Mandal, Katherine Klymko, and Michael R.DeWeese, “Entropy production and fluctuation theoremsfor active matter,” Phys. Rev. Lett. , 258001 (2017).[25] Patrick Pietzonka and Udo Seifert, “Entropy produc-tion of active particles and for particles in active baths,”J. Phys. A: Math. Theor. , 01LT01 (2017).[26] Cesare Nardini, ´Etienne Fodor, Elsen Tjhung, Fr´ed´ericvan Wijland, Julien Tailleur, and Michael E. Cates, “En-tropy production in field theories without time-reversalsymmetry: Quantifying the non-equilibrium character ofactive matter,” Phys. Rev. X , 021007 (2017).[27] Suraj Shankar and M. Cristina Marchetti, “Hidden en-tropy production and work fluctuations in an ideal activegas,” Phys. Rev. E , 020604 (2018).[28] Lennart Dabelow, Stefano Bo, and Ralf Eich-horn, “Irreversibility in active matter systems:Fluctuation theorem and mutual information,”Phys. Rev. X , 021009 (2019).[29] Shoichi Toyabe, Tetsuaki Okamoto, Takahiro Watanabe-Nakayama, Hiroshi Taketani, Seishi Kudo, and EiroMuneyuki, “Nonequilibrium energetics of a single f -atpase molecule,” Phys. Rev. Lett. , 198103 (2010).[30] Thomas Speck, “Stochastic thermodynamics for activematter,” EPL (Europhys. Lett.) , 30006 (2016).[31] ´E. Fodor, W. W. Ahmed, M. Almonacid, M. Bussonnier,N. S. Gov, M.-H. Verlhac, T. Betz, P. Visco, and F. vanWijland, “Nonequilibrium dissipation in living oocytes,”EPL (Europhys. Lett.) , 30008 (2016). [32] Laura Tociu, ´Etienne Fodor, Takahiro Nemoto,and Suriyanarayanan Vaikuntanathan, “How dissipa-tion constrains fluctuations in nonequilibrium liq-uids: Diffusion, structure, and biased interactions,”Phys. Rev. X , 041026 (2019).[33] ´Etienne Fodor, Takahiro Nemoto, and SuriyanarayananVaikuntanathan, “Dissipation controls transport andphase transitions in active fluids: mobility, diffusion andbiased ensembles,” New J. Phys. , 013052 (2020).[34] F. Cagnetta, F. Corberi, G. Gonnella, andA. Suma, “Large fluctuations and dynamic phasetransition in a system of self-propelled particles,”Phys. Rev. Lett. , 158002 (2017).[35] Trevor GrandPre and David T Limmer, “Currentfluctuations of interacting active brownian particles,”Phys. Rev. E , 060601 (2018).[36] Stephen Whitelam, Katherine Klymko, and DibyenduMandal, “Phase separation and large deviations of latticeactive matter,” J. Chem. Phys. , 154902 (2018).[37] Takahiro Nemoto, ´Etienne Fodor, Michael E. Cates,Robert L. Jack, and Julien Tailleur, “Optimizing ac-tive work: Dynamical phase transitions, collective mo-tion, and jamming,” Phys. Rev. E , 022605 (2019).[38] Giacomo Gradenigo and Satya N Majumdar, “Afirst-order dynamical transition in the displacementdistribution of a driven run-and-tumble particle,”J. Stat. Mech. , 053206 (2019).[39] Emil Mallmin, Richard A Blythe, and Martin R Evans,“A comparison of dynamical fluctuations of biased dif-fusion and run-and-tumble dynamics in one dimension,”J. Phys. A: Math. Theor. , 425002 (2019).[40] F. Cagnetta and E. Mallmin, “Efficiency of one-dimensional active transport conditioned on motility,”Phys. Rev. E , 022130 (2020).[41] Pietro Chiarantoni, Francesco Cagnetta, Federico Cor-beri, Giuseppe Gonnella, and Antonio Suma, “Work fluc-tuations of self-propelled particles in the phase separatedstate,” J. Phys. A: Math. Theor. (2020).[42] Trevor GrandPre, Katherine Klymko, Kranthi K. Man-dadapu, and David T. Limmer, “Entropy productionfluctuations encode collective behavior in active matter,”ArXiv e-prints (2020), arXiv:2007.12149.[43] V. Lecomte, C. Appert-Rolland, and F. van Wijland,“Thermodynamic formalism for systems with markov dy-namics,” J. Stat. Phys. , 51–106 (2007).[44] Hugo Touchette, “The large deviation approach to sta-tistical mechanics,” Phys. Rep. , 1–69 (2009).[45] Robert L. Jack and Peter Sollich, “Large deviationsand ensembles of trajectories in stochastic models,”Prog. Theor. Phys. Supp. , 304–317 (2010).[46] R. L. Jack and P. Sollich, “Effective interac-tions and large deviations in stochastic processes,”Eur. Phys. J. Special Topics , 2351–2367 (2015).[47] R. L. Jack, “Ergodicity and large deviationsin physical systems with stochastic dynamics,”Eur. Phys. J. B , 74 (2020).[48] Cristian Giardin`a, Jorge Kurchan, and Luca Peliti,“Direct evaluation of large-deviation functions,”Phys. Rev. Lett. , 120603 (2006).[49] Takahiro Nemoto, Freddy Bouchet, Robert L.Jack, and Vivien Lecomte, “Population-dynamicsmethod with a multicanonical feedback control,”Phys. Rev. E , 062123 (2016). [50] J. P. Garrahan, R. L. Jack, V. Lecomte, E. Pitard, K. vanDuijvendijk, and F. van Wijland, “Dynamical first-order phase transition in kinetically constrained modelsof glasses,” Phys. Rev. Lett. , 195702 (2007).[51] Lester O. Hedges, Robert L. Jack, Juan P. Garra-han, and David Chandler, “Dynamic order-disorderin atomistic models of structural glass formers,”Science , 1309–1313 (2009).[52] Thomas Speck, Alex Malins, and C. Patrick Roy-all, “First-order phase transition in a model glassformer: Coupling of local structure and dynamics,”Phys. Rev. Lett. , 195703 (2012).[53] Julien Tailleur and Jorge Kurchan, “Probing rare phys-ical trajectories with lyapunov weighted dynamics,”Nat. Phys. , 203 (2007).[54] Tanguy Laffargue, Khanh-Dang Nguyen ThuLam, Jorge Kurchan, and Julien Tailleur,“Large deviations of lyapunov exponents,”J. Phys. A: Math. Theor. , 254002 (2013).[55] Jakub Dolezal and Robert L Jack, “Large deviations andoptimal control forces for hard particles in one dimen-sion,” J. Stat. Mech. , 123208 (2019).[56] Salvatore Torquato and Frank H. Stillinger, “Local den-sity fluctuations, hyperuniformity, and order metrics,”Phys. Rev. E , 041113 (2003).[57] Raphal Chetrite and Hugo Touchette, “Variational andoptimal control representations of conditioned and drivenprocesses,” J. Stat. Mech. , P12001 (2015).[58] Ken Sekimoto, “Langevin equation and thermodynam-ics,” Prog. Theor. Phys. Supp. , 17–27 (1998).[59] Udo Seifert, “Stochastic thermodynamics, fluc-tuation theorems and molecular machines,”Rep. Prog. Phys. , 126001 (2012).[60] Andrea Puglisi and Umberto Marini Bettolo Marconi,“Clausius relation for active particles: What can we learnfrom fluctuations,” Entropy , 356 (2017).[61] Umberto Marini Bettolo Marconi, Andrea Puglisi, andClaudio Maggi, “Heat, temperature and clausius in-equality in a model for active brownian particles,”Sci. Rep. , 46496 (2017).[62] Vivien Lecomte and Julien Tailleur, “A numericalapproach to large deviations in continuous time,”J. Stat. Mech. , P03004–P03004 (2007).[63] Tobias Brewer, Stephen R Clark, Russell Brad-ford, and Robert L Jack, “Efficient characterisa-tion of large deviations using population dynamics,”J. Stat. Mech. , 053204 (2018).[64] AB Slowman, MR Evans, and RA Blythe, “Jammingand attraction of interacting run-and-tumble randomwalkers,” Phys. Rev. Lett. , 218101 (2016).[65] Rapha¨el Chetrite and Hugo Touchette, “Nonequilibriummarkov processes conditioned on large deviations,” inAnnales Henri Poincar´e, Vol. 16 (Springer, 2015) pp.2005–2057.[66] Daniel Jacobson and Stephen Whitelam, “Direct evalua-tion of dynamical large-deviation rate functions using avariational ansatz,” Phys. Rev. E , 052139 (2019).[67] Lorenzo Bertini, Alberto De Sole, DavideGabrielli, Giovanni Jona-Lasinio, and Clau-dio Landim, “Macroscopic fluctuation theory,”Rev. Mod. Phys. , 593–636 (2015).[68] M. E. Cates and J. Tailleur, “When are active brown- ian particles and run-and-tumble particles equivalent?consequences for motility-induced phase separation,”EPL (Europhys. Lett.) , 20010 (2013).[69] Mourtaza Kourbane-Houssene, Cl´ement Erignoux,Thierry Bodineau, and Julien Tailleur, “Exacthydrodynamic description of active lattice gases,”Phys. Rev. Lett. , 268003 (2018).[70] F. D. C. Farrell, M. C. Marchetti, D. Maren-duzzo, and J. Tailleur, “Pattern formation in self-propelled particles with density-dependent motility,”Phys. Rev. Lett. , 248101 (2012).[71] Robert L. Jack, Ian R. Thompson, and Peter Sol-lich, “Hyperuniformity and phase separation in bi-ased ensembles of trajectories for diffusive systems,”Phys. Rev. Lett. , 060601 (2015).[72] C. Appert-Rolland, B. Derrida, V. Lecomte, and F. vanWijland, “Universal cumulants of the current in diffusivesystems on a ring,” Phys. Rev. E , 021122 (2008).[73] S. Sachdev, Quantum Phase Transitions (CambridgeUniversity Press, 1999).[74] Alexandre P. Solon, Hugues Chat´e, and Julien Tailleur,“From phase to microphase separation in flocking mod-els: The essential role of nonequilibrium fluctuations,”Phys. Rev. Lett. , 068101 (2015).[75] Yann-Edwin Keta, “yketa/active work,” GitHub reposi-tory (2019), MIT licensed.[76] Hugo Touchette, “Introduction to dynamical large devia-tions of markov processes,” Physica A , 5–19 (2018).[77] Lars Onsager and Stefan Machlup, “Fluctuations and ir-reversible processes,” Phys. Rev. , 1505 (1953).[78] Leticia F Cugliandolo and Vivien Lecomte, “Rules of cal-culus in the path integral representation of white noiselangevin equations: the onsager–machlup approach,”J. Phys. A: Math. Theor. , 345001 (2017).[79] Leticia F Cugliandolo, Vivien Lecomte, andFr´ed´eric van Wijland, “Building a path-integralcalculus: a covariant discretization approach,”J. Phys. A: Math. Theor. , 50LT01 (2019).[80] Avishek Das and David T. Limmer, “Varia-tional control forces for enhanced sampling ofnonequilibrium molecular dynamics simulations,”J. Chem. Phys. , 244123 (2019).[81] Ushnish Ray, Garnet Kin-Lic Chan, and David T.Limmer, “Exact fluctuations of nonequilibriumsteady states from approximate auxiliary dynamics,”Phys. Rev. Lett. , 210602 (2018).[82] Thibaut Arnoulx de Pirey, Gustavo Lozano, and Fr´ed´ericvan Wijland, “Active hard spheres in infinitely many di-mensions,” Phys. Rev. Lett. , 260602 (2019).[83] Arghya Das, Abhishek Dhar, and Anupam Kundu, “Gapstatistics of two interacting run and tumble particles inone dimension,” J. Phys. A: Math. Theor. (2020).[84] AB Slowman, MR Evans, and RA Blythe, “Ex-act solution of two interacting run-and-tumblerandom walkers with finite tumble duration,”J. Phys. A: Math. Theor.50