# Out of equilibrium response and fluctuation-dissipation violations across
scales in flocking systems

Federica Ferretti  Current address: NSF-Simons National Institute for Theory
and Mathematics in Biology, 60611 Chicago, IL  Department of Chemical
Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts
02139, USA  Irene Giardina<EMAIL_ADDRESS>Dipartimento di Fisica,
Universit‘a Sapienza, 00185 Rome, Italy  Istituto Sistemi Complessi, Consiglio
Nazionale delle Ricerche, UOS Sapienza, 00185 Rome, Italy  INFN, Unit‘a di
Roma 1, 00185 Rome, Italy  Tomas Grigera  Instituto de Física de Líquidos y
Sistemas Biológicos (IFLySiB), Universidad Nacional de La Plata and CONICET,
La Plata, Argentina  CCT CONICET La Plata, Consejo Nacional de Investigaciones
Científicas y Técnicas, La Plata, Argentina  Departamento de Física, Facultad
de Ciencias Exactas, Universidad Nacional de La Plata, La Plata, Argentina
Istituto Sistemi Complessi, Consiglio Nazionale delle Ricerche, UOS Sapienza,
00185 Rome, Italy  Giulia Pisegna  Max Planck Institute for Dynamics and Self-
Organization (MPI-DS), D-37077 Goettingen, Germany  Mario Veca
<EMAIL_ADDRESS>Dipartimento di Fisica, Universit‘a Sapienza, 00185
Rome, Italy

###### Abstract

Flocking systems are known to be strongly out of equilibrium. Energy input
occurs at the individual level to ensure self-propulsion, and the individual
motility in turn contributes to ordering, enhancing information propagation
and strengthening collective motion. However, even beyond ordering, a crucial
feature of natural aggregations is response. How, then, do off-equilibrium
features affect the response of the system? In this work, we consider a
minimal model of flocking and investigate response behavior under directional
perturbations. We show that equilibrium dynamical fluctuation-dissipation
relations between response and correlations are violated, both at the local
and at the global level. The amount of violation peaks at the ordering
transition, exactly as for the entropy production rate. Entropy is always
produced locally and connected to the local fluctuation-dissipation violation
via Harada-Sasa relationships. However, cooperative mechanisms close to the
transition spread off-equilibrium effects to the whole system, producing an
out of equilibrium response on the global scale. Our findings elucidate the
role of activity and interactions in the cost repartition of collective
behavior and explain what observed in experiments on natural living groups.

Introduction. Collective motion in living systems has been observed on
different scales and regimes, from swimming bacteria and cell colonies to
insect swarms and bird flocks Krause and Ruxton (2002). From the perspective
of statistical physics, these phenomena represent paradigmatic instances of
non-equilibrium emergent behavior and they have been intensively studied, both
at theoretical and experimental level, in the context of active matter
Marchetti  _et al._ (2013); Ramaswamy (2010); Gompper  _et al._ (2020).
Broadly speaking, a flock has the following distinctive features: individuals
are ‘active’, i.e. they are endowed with a self-propulsion mechanism
transforming energy into motion, and they coordinate with each other via short
range interaction rules. Individual motility is responsible for the non
equilibrium character of the problem. Interactions produce collective patterns
on the large scale. The combination of these two ingredients leads to a non-
trivial phenomenology comprising kinetic ordering transitions, novel classes
of critical behavior and motility induced phase separation Vicsek  _et al._
(1995); Vicsek and Zafeiris (2012); Chaté  _et al._ (2008); Cates and Tailleur
(2015).

The impact of activity on order has been investigated at length. Seminal works
Vicsek  _et al._ (1995); Toner and Tu (1995, 1998) have indeed shown that
self-propulsion favours ordering, by enhancing information propagation through
the system and stabilizing the transition to collective motion. On the other
hand, much less is known on how non-equilibrium features affect the response
of the system to perturbations. Response is a key trait for survival in living
aggregations, and it therefore represents a quantity of primary biological
relevance. Besides, even more than order itself, it provides a bona fide
signature of collective behavior, by quantifying the ability to retain global
coherence in presence of external signals or threats. In this work, we address
systematically this question. We consider an archetypal model of flocking and
compute analytically and numerically the response function. We show that
equilibrium dynamical fluctuation-dissipation relations between response and
correlations are violated on all timescales, and we quantify the corresponding
deviation. The global amount of the non-equilibrium contribution depends on
the noise and it peaks at the transition point. We explain how non-equilibrium
effects distribute across scales: energy is injected at the individual level
via self-propulsion and violation of detailed balance occurs on the scale of
the interaction range; however, close to the transition cooperative mechanisms
arise in the dynamics and the non-equilibrium effects spread on all scales.

The model and the perturbation protocol. We consider a continuous time
generalization of the Vicsek model of flocking Vicsek  _et al._ (1995);
Chepizhko  _et al._ (2021) where self-propelled particles regulate their own
speed using a speed control potential, and interact with each other locally
via mutual adjustment of their velocities Bialek  _et al._ (2014); Cavagna
_et al._ (2022). The model captures in minimal terms the imitative nature of
behavioral coordination rules in aggregations of living organisms Krause and
Ruxton (2002), as well as the effective interactions present in many non-
living active assemblies Vicsek and Zafeiris (2012). The equations of motion
are

$\frac{d\mathbf{v}_{i}}{dt}=-\frac{\partial\mathcal{H}}{\partial\mathbf{v}_{i}}+\bm{\xi}_{i}\
\ \ ,\ \ \frac{d\mathbf{r}_{i}}{dt}=\mathbf{v}_{i}$ (1)

where $\bm{\xi}_{i}$ is a white noise with variance
$\langle\xi_{i\alpha}(t)\xi_{j\beta}(t^{\prime})\rangle=2T\delta_{ij}\delta_{\alpha\beta}\delta(t-t^{\prime})$,
and the pseudo Hamiltonian ${\mathcal{H}}$ reads

$\mathcal{H}=\frac{J}{2}\sum_{ij}n_{ij}(\mathbf{v}_{i}-\mathbf{v}_{j})^{2}+\frac{g}{2}\sum_{i}(|\mathbf{v}_{i}|-v_{0})^{2}\
.$ (2)

The first term in this expression describes mutual adjustment between velocity
vectors, $J$ being the scale of the interaction and $n_{ij}(t)$ the adjacency
matrix. In the following, we will adopt a metric interaction rule, thereby
$n_{ij}(t)=1$ if particles $i$ and $j$ at time $t$ have mutual distance
smaller than a fixed range $r_{0}=1$, and $n_{ij}(t)=0$ otherwise. The second
term in Eq. (2) is the speed control, setting the individual speeds close to a
cruising reference value $v_{0}$, with fluctuations regulated by the parameter
$g$ 111 We could have chosen different forms for the speed potential (e.g.
functions of ${\mathbf{v}}_{i}^{2}$ rather than $|\mathbf{v}_{i}|$). In the
region around the ordering transition - which is the most interesting regime
for our analysis - there is no much difference, so we stick to the simple
Gaussian shape of Eq. 2. . The fixed modulus Vicsek model is recovered for
$g\to\infty$.

With decreasing noise strength, model (1) exhibits a kinetic ordering
transition to a state of collective motion, where the mean group velocity
${\mathbf{V}}=(1/N)\sum_{i}{\mathbf{v}}_{i}$ and the global degree of
alignment - the polarization
$\bm{\Phi}=(1/N)\sum_{i}\mathbf{v}_{i}/|\mathbf{v}_{i}|$ \- are different from
zero (see App. A). Our aim here is to investigate response behavior, and to
understand how it relates to ordering and to the out-of-equilibrium features
of the model. To do so, we consider a perturbation protocol where external
fields are applied to the individual velocities, i.e.
$\mathcal{H}\to\mathcal{H}-\sum_{i}{\mathbf{h}}_{i}\cdot{\mathbf{v}}_{i}$. The
same protocol has been recently used in Chaté  _et al._ (2008); Kyriakopoulos
_et al._ (2016); Brambati  _et al._ (2022); Loffredo  _et al._ (2023));
however, no systematic investigation has been performed so far of fluctuation-
dissipation relations.

Response functions and fluctuation-dissipation relations. The linear response
functions can be defined in the usual way as
$R^{\alpha\beta}_{ij}(t,s)=\delta\langle v_{i}^{\alpha}(t)\rangle_{h}/\delta
h_{j}^{\beta}(s)|_{h=0}$ Marconi  _et al._ (2008). In equilibrium systems,
such response functions are related to the velocity correlation functions
$C^{\alpha\beta}_{ij}(t,s)=\langle v_{i}^{\alpha}(t)v_{j}^{\beta}(s)\rangle$
by the dynamical Fluctuation Dissipation Theorem (FDT) Kubo (1966). However,
model (1) is intrinsically out of equilibrium due to the active nature of the
particles, and we expect violations of such relationships. In the stationary
state, for the diagonal spacial components (the only ones that, for symmetry
reasons, are non-zero in our case), we can then write (for $t>s)$

$TR^{\alpha\alpha}_{ij}(t-s)=-\frac{dC^{\alpha\alpha}_{ij}(t-s)}{dt}-d^{\alpha\alpha}_{ij}(t-s)$
(3)

where $d^{\alpha\alpha}_{ij}$ measures the deviation from equilibrium.

In general, computing either experimentally or numerically the full response
tensor is a demanding task, as one would need to apply local fields and
monitor reactions for all pairs of individuals. On the contrary, global
responses are more easily accessible and statistically more reliable. Let us
then consider the simpler setup where we apply a uniform field
$h_{i}^{\alpha}(t)=h_{\alpha}(t)$ on the system and we measure how the global
mean velocity ${\mathbf{V}}$ is affected by the perturbation. In this case,
the global linear response is given by $R_{\alpha\alpha}(t-s)=\delta\langle
V_{\alpha}(t)\rangle/\delta
h_{\alpha}(s)|_{h=0}=(1/N)\sum_{ij}R^{\alpha\alpha}_{ij}(t-s)$, and the
analogue of Eq. (3) is

$TR_{\alpha\alpha}(t-s)=-\frac{dC_{\alpha\alpha}(t-s)}{dt}-d_{\alpha\alpha}(t-s)\
,$ (4)

with $C_{\alpha\alpha}(t-s)=N\langle V_{\alpha}(t)V_{\alpha}(s)\rangle$ and
$d_{\alpha\alpha}(t-s)=(1/N)\sum_{ij}d_{ij}^{\alpha\alpha}(t-s)$.

Interestingly, for Langevin dynamics the response functions can be explicitly
computed starting from the equations of motion. To this task, it is convenient
to introduce the Onsager-Machlup action Onsager and Machlup (1953)
$S[\\{\mathbf{v}_{i},\mathbf{r}_{i}\\}]$, defined by

$\displaystyle
P\left[\\{\mathbf{v}_{i}(t),\bf{r}_{i}(t)\\}\right]=\frac{1}{\cal
N}\exp{\\{-S[\\{\mathbf{v}_{i}(t),\mathbf{r}_{i}(t)\\}]\\}}$ (5)
$\displaystyle\quad=\biggl{\langle}\delta\left({\mathbf{v}}_{i}(t)-{\mathbf{v}}_{i}(0)-\int_{0}^{t}dt^{\prime}\
[-\partial{\cal
H}/\partial{\mathbf{v}}_{i}+{\bm{\xi}}_{i}(t^{\prime})]\right)\biggr{\rangle}\
.$

where averages are performed over the trajectory realizations. In terms of the
action $S$, the response functions can be written as
$R^{\alpha\alpha}_{ij}(t,s)=-\langle v_{i}^{\alpha}(t)\ \delta S/\delta
h_{j}^{\alpha}(s)\rangle|_{h=0}$. Using standard techniques of stochastic
calculus Cugliandolo (2011); Marconi  _et al._ (2008); Caprini (2021) one then
finds, for $t>s$, (see App. B for details)

$d_{\alpha\alpha}=\frac{gNv_{0}}{2}[\langle
V_{\alpha}(t)\Phi_{\alpha}(s)\rangle-\langle
V_{\alpha}(s)\Phi_{\alpha}(t)\rangle]\ .$ (6)

The deviation from equilibrium is now expressed in terms of the global
correlation between velocity and polarization, computed in zero field. Eq. (6)
therefore allows to precisely quantify the FDT violation and to compute the
response (via Eq. (4)) without actually perturbing the system. We note that in
the r.h.s of Eq. (6) the speed control parameter $g$ explicitly appears. In
fact, as we will show later on, the deviation only mildly depends on $g$ and a
well defined asymptotic $g\to\infty$ limit is quickly reached as soon as
$g\sim{\rm O}(1)$. Still, allowing the individual speed to fluctuate is
crucial to derive Eq. (6). In other terms, speed fluctuations are not
indispensable to produce out-of-equilibrium effects, but they can be used very
effectively to capture and quantify them (see App.E for more discussion).

FDT violations across the phase diagram. We now proceed to compute numerically
the response functions and the deviations from equilibrium, and monitor them
across the phase diagram of the model. We use a discretized version of Eqs.
(1), in $d=3$, following a standard Euler integration scheme. We select a
value of density for which violations are clearly visible in the whole phase
space, and then tune the noise strength to span from the disordered to the
ordered region. We focus initially on checking the validity of Eq.(6). To do
so, we compute the response in two ways, either perturbing the system or using
Eqs. (4) and (6) in absence of perturbation. In fact, it is numerically more
convenient to look at the integrated version of Eq. (4), i.e.
$T\chi_{\alpha\alpha}(t)=[C_{\alpha\alpha}(0)-C_{\alpha\alpha}(t)]-D_{\alpha\alpha}(t)$,
where $\chi_{\alpha\alpha}(t)=\int_{0}^{t}d\tau R_{\alpha\alpha}(\tau)$ is the
integrated response Cugliandolo (2011) and
$D_{\alpha\alpha}(t)=\int_{0}^{t}d\tau d_{\alpha\alpha}(\tau)$. For $T>T_{c}$
all coordinates $\alpha$ are statistically equivalent; for $T<T_{c}$ we
consider the longitudinal component, i.e. we choose $\alpha$ as the direction
of the polarization (see App. A for details). The result is displayed in the
inset of Fig.1b, where we can clearly see that the response curves obtained in
the two ways coincide within statistical fluctuations. Once convinced that we
can trust Eq.(6), we can then exploit this expression to compute the deviation
function, and therefore the response, without actually perturbing the system.
The deviation function $d(t)$ is displayed in Fig.1a for several values of the
temperature (we drop the coordinate index for notation convenience). In
general, $d(t)$ shows a first quick increase on short timescales towards a
maximum, and a slower decay to zero at longer times. We shall comment later on
the origins of such dynamical behavior. For the time being, we notice that the
violation is larger and more extended in time as the temperature approaches
the critical value $T_{c}$ of the kinetic transition. A different way - which
is standard in the literature - to capture out-of-equilibrium behavior is to
plot the integrated response parametrically as a function of the correlation
function Bouchaud  _et al._ (1998). In equilibrium dynamics the expected
behaviour is a straight line with slope given by (minus) the inverse
temperature. On the contrary, as shown in Fig.1b, in our model the curves
clearly depart from the equilibrium one, and such departure looks more
pronounced close to criticality, consistently with Fig.1a.

Figure 1: FDT Violation. a) Deviation $d(t)$ as a function of time, for
different values of the temperature. b) Integrated response vs. correlation,
for different values of the temperature; both quantities are normalized by
$C(t=0)$. Inset: integrated response computed by perturbing the system (red
points) or by using Eqs. (6) and (4) (blue points), for $T=15$. c) Maximum of
the deviation function, $d_{\rm max}$, as a function of temperature. Inset:
behavior of $d_{\rm max}$ as a function of $g$. d) Asymptotic off-equilibrium
offset $T_{\rm eff}(\infty)/T-1$ as a function of temperature; the black line
corresponds to equilibrium. The parameters of the simulation are
$J=1.5,v_{0}=4,g=4,\rho=0.25,N=512$, with temperatures described by the color
bar and $T_{c}=10$ (light green).

To better investigate the role of temperature, we now introduce two synthetic
quantifiers of the departure from equilibrium. As a first measure, we consider
the maximum $d_{\rm max}$ of the deviation function $d(t)$, which occurs on
short timescales. In Fig. 1c $d_{\rm max}$ is plotted as a function of
temperature, and it displays a sharp maximum at the transition temperature. In
the inset we also show $d_{\rm max}(T_{c})$ as a function of $g$: as
anticipated, it quickly converges to the limit $g\to\infty$ (i.e. standard
Vicsek).

Another quantity that is widely used in the literature to quantify FDT
violations is the so-called effective temperature Hohenberg and Shraiman
(1989); Cugliandolo  _et al._ (1997); Cugliandolo (2011), defined as the ratio
between the correlation and the integrated response $T_{\rm
eff}(t)=(C(0)-C(t))/\chi(t)$. Even though $T_{\rm eff}$ represents a real
temperature only in a few cases, we will stick for convenience to this
terminology. In Fig. 1d we display the asymptotic value $(T_{\rm
eff}(t\to\infty)-T)/T$ (which is zero at equilibrium) as a function of
temperature. Also this quantity - which encodes FDT violations on long
timescales - displays a maximum in correspondence of the critical temperature.

FDT violations across scales. The above results indicate that FDT violations
are maximal at the transition, and that the response is strongly out of
equilibrium on all time-scales. The behavior in time of the deviation function
$d(t)$, however, suggests that different processes might regulate the short
time regime and the long time decay. In Fig. 2a we plot the self contribution
$d_{\rm self}=(1/N)\sum_{i}d_{ii}$ to the deviation and compare it to the full
function. This contribution dominates the first regime of $d(t)$ and it decays
to zero on a timescale corresponding to the location of $d_{\rm max}$.

Interestingly, $d_{\rm self}$ is connected to another quantifier of non-
equilibrium, the Entropy Production Rate (EPR), which measures the amount of
violation of detailed balance in the system. The entropy production is defined
as the logarithm of the ratio between the probability of a forward and a
backward trajectory Seifert (2012); Ritort (2008), and it can be computed from
the Onsager-Machlup action. We get

${\rm
EPR}=\lim_{\tau\to\infty}\frac{1}{\tau}\Big{\langle}\frac{J}{2T}\int_{0}^{\tau}d\tau\sum_{ij}\frac{dn_{ij}}{dt}\circ\sum_{\alpha}(v^{\alpha}_{i}-v^{\alpha}_{j})^{2}\Big{\rangle}$
(7)

where $\circ$ denotes the Stratonovitch prescription (see App.C). A
discretized version of this expression can be used to evaluate the EPR from
numerical trajectories. The result is displayed in Fig.2b, and it shows that
the EPR also peaks at the critical temperature Ferretti  _et al._ (2022); Yu
and Tu (2022). In fact, the EPR and FDT violations are directly related to
each other via Harada-Sasa relationships Harada and Sasa (2005), which we
explicitly derive in App. D,

$\displaystyle{\rm EPR}$
$\displaystyle=\sum_{i,\alpha}\int\frac{d\omega}{2\pi}\frac{\omega}{T}\left[\omega\tilde{C}_{ii}^{\alpha\alpha}(\omega)-2T\mathfrak{Im}(\tilde{R}_{ii}^{\alpha\alpha}(\omega))\right]$
$\displaystyle=\frac{N}{T}\frac{d}{dt}\Big{(}\sum_{\alpha}d_{\rm
self}^{\alpha\alpha}(t)\Big{)}\Big{|}_{t=0^{+}}$ (8)

We recognize in the brackets the Fourier transform of the elements
$d^{\alpha\alpha}_{ii}(t)$ of the deviation matrix, which also contribute to
$d_{\rm self}$. We have verified this relation, by computing Eq. (8) via
Mallivian weight sampling Warren and Allen (2012) and comparing the result
with Eq. (7), see Fig.2b and App. D.

The connection with the EPR helps understanding the origin of the short time
regime of the FDT violation: it is directly determined by the production of
entropy in the system. Eq.(7) shows that the EPR depends on the rate of change
of the interaction network $dn_{ij}/dt$, that is on how quickly nearby
particles enter and exit their interaction range. It is therefore a local
process. Besides, since mutual diffusion in space is maximal at the transition
Ferretti  _et al._ (2022); Tu  _et al._ (1998); Chaté  _et al._ (2008) the EPR
peaks at $T=T_{c}$ 222At least for not too small values of $g$. We remind that
$g$ regulates the amplitude of speed fluctuations. If it is too small speeds
are not confined and they can alter the scenario so far described. This is
however an unrealistic regime, which is not of interest to our analysis..

As we already noticed, strong deviations from equilibrium are observed in the
response also on the long timescales, as testified by the slow decay of the
function $d(t)$ and by the behavior of $T_{\rm eff}(t=\infty)$ across the
phase diagram (see Fig. 1d). To explain this behavior, it is crucial to
consider the interacting nature of the system and, in particular, the
occurrence of strong cooperative phenomena as the transition point is
approached. At large times, the deviation function is dominated by the non-
self contribution, and it therefore involves correlations between pairs of
particles. Close to $T_{c}$ such correlations become long-range, distributing
off-equilibrium effects to the large scales and slowing down the decay of
$d(t)$.

To confirm the collective nature of FDT violations at large times, we
performed simulations for systems of different size $N$, at the critical
temperature $T=T_{c}(N)$. The resulting $d(t)$ curves are displayed in Fig.
2c: the short time increase of the deviation function, being related to a
local process, is always the same; on the contrary, the long time decay is
slower the larger the size, a behavior typical of critical slowing down.
Indeed, we show in Fig. 2d that the deviation curves obey dynamic scaling,
falling one on top of the other once the time is rescaled by a characteristic
time growing with the system’s size as $L^{z}$. The value of the dynamic
exponent used for the collapse is $z=1.7$, that is the value predicted by the
incompressible field theory of the Vicsek model Chen  _et al._ (2015). We
conclude that - as long as the system is homogeneous enough in density and the
incompressible regime holds Di Carlo and Scandolo (2022) \- long time FDT
violations are a critical effect. Above a (density dependent) crossover size
the cooperative processes at the transition are dominated by heterogeneities
and phase separation Chaté  _et al._ (2008); Di Carlo and Scandolo (2022) and
we expect the long time dynamics to be ruled by different laws. The
simulations performed in this paper, as well as experiments on natural groups
in $d=3$ Cavagna  _et al._ (2008, 2023), all belong to the incompressible
regime (see also App.A.).

Figure 2: FDT violations across scales. a) Self-contribution $d_{\rm self}$
(red points) and full deviation (blue points) as a function of time, at
$T=T_{c}$. b) EPR computed using the trajectories with Eq. (7) (red points)
and the Harada-Sasa formula (blue points). c) Deviation curves $d(t)$ at
$T=T_{c}(N)$ for several values of the system’s size. d) Rescaled deviation
$d(t)/d_{\rm max}(N)$ as a function of the rescaled time $tL^{-z}$, with
$z=1.7$. The parameters of the simulation are
$J=1.5,v_{0}=4,g=4,\rho=0.25,N=512$ (a and b). Results at $\rho=1.5$ can be
found in App. A .

The general picture can be summarized as follows. Energy is injected in the
system at the individual level, via the particle self-propulsion and the
relative speed control function. Then, entropy is produced locally, every time
two particles start/stop interacting with each other, and it is maximal at
$T_{c}$ due to the fast network reshuffling. Violations of detailed balance
immediately alter FDT relations on the short scales. The existence of long-
range cooperative processes close to the transition then propagate off-
equilibrium effects through the whole system. FDT deviations therefore persist
at large times, and their overall amount peaks at the transition point.

Conclusions. The way a system manages its energy budget is a crucial issue for
living systems, which typically acquire external resources to build functional
ordered structures and dissipate excess entropy in the environment.
Understanding how this occurs and what determines the repartition of this
process is therefore an important question Gnesotto  _et al._ (2018); Lynn
_et al._ (2022). Indeed, it remains in general unclear whether activity in the
individuals of an aggregation affects the dynamics also at the macroscopic
level Egolf (2000); for example, recent analysis indicate that in motility
induced phase separation equilibrium is restored at large scales Nardini  _et
al._ (2017); Maggi  _et al._ (2022). Our results show that in fact in flocking
systems non-equilibrium features do manifest on all scales, being stronger on
the short times, and slowly decaying at long ones. The key ingredient for this
to occur is the combination of the individual driving enforcing self-
propulsion, and the presence of an ordering transition coupling dynamical
modes over the whole system.

Our analysis indicates that FDT violations give a more complete description of
out-of equilibrium behavior than the EPR, as they capture violations on all
scales, whereas the EPR is a strictly a local quantity. The different role of
EPR and FDT violations has been pointed out in several works both at
microscopic and at field theory level Crisanti  _et al._ (2012); Gnesotto  _et
al._ (2018); Nardini  _et al._ (2017); Caballero and Cates (2020), but it
becomes prominent in strongly interacting systems as the one considered here.
In some cases, it might be useful to define a coarse-grained EPR Esposito
(2012) capturing dissipation at different scales. In this way, an inverse
dissipation scaling similar to what described here has been recently observed
in self-similar reaction networks Yu  _et al._ (2021), self-similarity in the
network structure playing the role of the scale-free behavior at the critical
point.

Finally, we note that our results explain what is observed in experiments on
natural flocks and swarms: while flocks of birds behave as quasi-equilibrium
systems Mora  _et al._ (2016), swarms of midges exhibit strongly out of
equilibrium features Cavagna  _et al._ (2023). This is consistent with our
findings: flocks are highly polarized groups and therefore belong to a region
where both EPR and FDT violations are small (see Fig. 1); swarms, on the
contrary, are quasi-critical systems and obey dynamic scaling Attanasi  _et
al._ (2014); Cavagna  _et al._ (2017) thereby representing instances at the
peak of the violation-noise diagram.

###### Acknowledgements.

This work was supported by MIUR grant PRIN-2020PFCXPE and by ERC grant RG.BIO
(n. 785932). We thank A. Cavagna and J. Cristín for many interesting
discussions.

## References

  * Krause and Ruxton (2002) J. Krause and G. D. Ruxton,  _Living in groups_ (Oxford University Press, 2002).
  * Marchetti  _et al._ (2013) M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Rev. Mod. Phys. 85, 1143 (2013).
  * Ramaswamy (2010) S. Ramaswamy, Annual Review of Condensed Matter Physics 1, 323 (2010).
  * Gompper  _et al._ (2020) G. Gompper, R. G. Winkler, T. Speck, A. Solon, C. Nardini, F. Peruani, H. Löwen, R. Golestanian, U. B. Kaupp, L. Alvarez, _et al._ , Journal of Physics: Condensed Matter 32, 193001 (2020).
  * Vicsek  _et al._ (1995) T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Phys. Rev. Lett. 75, 1226 (1995).
  * Vicsek and Zafeiris (2012) T. Vicsek and A. Zafeiris, Physics Reports 517, 71 (2012).
  * Chaté  _et al._ (2008) H. Chaté, F. Ginelli, G. Grégoire, and F. Raynaud, Phys. Rev. E 77, 046113 (2008).
  * Cates and Tailleur (2015) M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys. 6, 219 (2015).
  * Toner and Tu (1995) J. Toner and Y. Tu, Phys. Rev. Lett. 75, 4326 (1995).
  * Toner and Tu (1998) J. Toner and Y. Tu, Phys. Rev. E 58, 4828 (1998).
  * Chepizhko  _et al._ (2021) O. Chepizhko, D. Saintillan, and F. Peruani, Soft Matter 17, 3113 (2021).
  * Bialek  _et al._ (2014) W. Bialek, A. Cavagna, I. Giardina, T. Mora, O. Pohl, E. Silvestri, M. Viale, and A. M. Walczak, Proceedings of the National Academy of Sciences 111, 7212 (2014).
  * Cavagna  _et al._ (2022) A. Cavagna, A. Culla, X. Feng, I. Giardina, T. S. Grigera, W. Kion-Crosby, S. Melillo, G. Pisegna, L. Postiglione, and P. Villegas, Nature Communications 13, 2315 (2022).
  * Note (1) We could have chosen different forms for the speed potential (e.g. functions of ${\mathbf{v}}_{i}^{2}$ rather than $|\mathbf{v}_{i}|$). In the region around the ordering transition - which is the most interesting regime for our analysis - there is no much difference, so we stick to the simple Gaussian shape of Eq. 2.
  * Kyriakopoulos  _et al._ (2016) N. Kyriakopoulos, F. Ginelli, and J. Toner, New Journal of Physics 18, 073039 (2016).
  * Brambati  _et al._ (2022) M. Brambati, G. Fava, and F. Ginelli, Physical Review E 106, 024608 (2022).
  * Loffredo  _et al._ (2023) E. Loffredo, D. Venturelli, and I. Giardina, Physical Biology 20, 035003 (2023).
  * Marconi  _et al._ (2008) U. M. B. Marconi, A. Puglisi, L. Rondoni, and A. Vulpiani, Physics reports 461, 111 (2008).
  * Kubo (1966) R. Kubo, Reports on progress in physics 29, 255 (1966).
  * Onsager and Machlup (1953) L. Onsager and S. Machlup, Physical Review 91, 1505 (1953).
  * Cugliandolo (2011) L. F. Cugliandolo, Journal of Physics A: Mathematical and Theoretical 44, 483001 (2011).
  * Caprini (2021) L. Caprini, Journal of Statistical Mechanics: Theory and Experiment 2021, 063202 (2021).
  * Bouchaud  _et al._ (1998) J.-P. Bouchaud, L. Cugliandolo, J. Kurchan, and M. Mézard, Directions in Condensed Matter Physics 12, 443 (1998).
  * Hohenberg and Shraiman (1989) P. Hohenberg and B. I. Shraiman, Physica D: Nonlinear Phenomena 37, 109 (1989).
  * Cugliandolo  _et al._ (1997) L. F. Cugliandolo, J. Kurchan, and L. Peliti, Physical Review E 55, 3898 (1997).
  * Seifert (2012) U. Seifert, Reports on progress in physics 75, 126001 (2012).
  * Ritort (2008) F. Ritort, (2008).
  * Ferretti  _et al._ (2022) F. Ferretti, S. Grosse-Holz, C. Holmes, J. L. Shivers, I. Giardina, T. Mora, and A. M. Walczak, Physical Review E 106, 034608 (2022).
  * Yu and Tu (2022) Q. Yu and Y. Tu, Physical review letters 129, 278001 (2022).
  * Harada and Sasa (2005) T. Harada and S.-i. Sasa, Physical review letters 95, 130602 (2005).
  * Warren and Allen (2012) P. B. Warren and R. J. Allen, Physical review letters 109, 250601 (2012).
  * Tu  _et al._ (1998) Y. Tu, J. Toner, and M. Ulm, Phys. Rev. Lett. 80, 4819 (1998).
  * Chaté  _et al._ (2008) H. Chaté, F. Ginelli, G. Grégoire, F. Peruani, and F. Raynaud, The European Physical Journal B 64, 451 (2008).
  * Note (2) At least for not too small values of $g$. We remind that $g$ regulates the amplitude of speed fluctuations. If it is too small speeds are not confined and they can alter the scenario so far described. This is however an unrealistic regime, which is not of interest to our analysis.
  * Chen  _et al._ (2015) L. Chen, J. Toner, and C. F. Lee, New Journal of Physics 17, 042002 (2015).
  * Di Carlo and Scandolo (2022) L. Di Carlo and M. Scandolo, New Journal of Physics (2022).
  * Cavagna  _et al._ (2008) A. Cavagna, A. Cimarelli, I. Giardina, A. Orlandi, G. Parisi, A. Procaccini, R. Santagati, and F. Stefanini, Mathematical biosciences 214, 32 (2008).
  * Cavagna  _et al._ (2023) A. Cavagna, L. Di Carlo, I. Giardina, T. S. Grigera, S. Melillo, L. Parisi, G. Pisegna, and M. Scandolo, Nature Physics 19, 1043 (2023).
  * Gnesotto  _et al._ (2018) F. S. Gnesotto, F. Mura, J. Gladrow, and C. P. Broedersz, Reports on Progress in Physics 81, 066601 (2018).
  * Lynn  _et al._ (2022) C. W. Lynn, C. M. Holmes, W. Bialek, and D. J. Schwab, Physical review letters 129, 118101 (2022).
  * Egolf (2000) D. A. Egolf, Science 287, 101 (2000).
  * Nardini  _et al._ (2017) C. Nardini, É. Fodor, E. Tjhung, F. Van Wijland, J. Tailleur, and M. E. Cates, Physical Review X 7, 021007 (2017).
  * Maggi  _et al._ (2022) C. Maggi, N. Gnan, M. Paoluzzi, E. Zaccarelli, and A. Crisanti, Communications Physics 5, 55 (2022).
  * Crisanti  _et al._ (2012) A. Crisanti, A. Puglisi, and D. Villamaina, Physical Review E 85, 061127 (2012).
  * Caballero and Cates (2020) F. Caballero and M. E. Cates, Phys. Rev. Lett. 124, 240604 (2020).
  * Esposito (2012) M. Esposito, Physical Review E 85, 041125 (2012).
  * Yu  _et al._ (2021) Q. Yu, D. Zhang, and Y. Tu, Physical review letters 126, 080601 (2021).
  * Mora  _et al._ (2016) T. Mora, A. M. Walczak, L. Del Castello, F. Ginelli, S. Melillo, L. Parisi, M. Viale, A. Cavagna, and I. Giardina, Nature Physics 12, 1153 (2016).
  * Attanasi  _et al._ (2014) A. Attanasi, A. Cavagna, L. Del Castello, I. Giardina, S. Melillo, L. Parisi, O. Pohl, B. Rossaro, E. Shen, E. Silvestri, and M. Viale, Phys. Rev. Lett. 113, 238102 (2014).
  * Cavagna  _et al._ (2017) A. Cavagna, D. Conti, C. Creato, L. Del Castello, I. Giardina, T. S. Grigera, S. Melillo, L. Parisi, and M. Viale, Nature Physics 13, 914 (2017).
  * Brézin  _et al._ (1973) E. Brézin, D. Wallace, and K. G. Wilson, Physical Review B 7, 232 (1973).
  * Brézin and Wallace (1973) E. Brézin and D. Wallace, Physical Review B 7, 1967 (1973).
  * Fischer  _et al._ (2018) L. P. Fischer, P. Pietzonka, and U. Seifert, Physical Review E 97, 022143 (2018).
  * Lebowitz and Spohn (1999) J. L. Lebowitz and H. Spohn, Journal of Statistical Physics 95, 333 (1999).
  * Warren and Allen (2013) P. B. Warren and R. J. Allen, Entropy 16, 221 (2013).
  * Cavagna  _et al._ (2014) A. Cavagna, I. Giardina, F. Ginelli, T. Mora, D. Piovani, R. Tavarone, and A. M. Walczak, Physical Review E 89, 042707 (2014).
  * Speck and Seifert (2006) T. Speck and U. Seifert, Europhysics Letters 74, 391 (2006).
  * Chetrite  _et al._ (2008) R. Chetrite, G. Falkovich, and K. Gawedzki, Journal of Statistical Mechanics: Theory and Experiment 2008, P08005 (2008).

##  APPENDIX A Numerical simulations

The equations of motions for our model are,

$\dfrac{d{\mathbf{r}}_{i}}{dt}=\mathbf{v}_{i},\quad\quad\quad\quad\dfrac{d{\mathbf{v}}_{i}}{dt}=-J\sum_{j}n_{ij}(t)(\mathbf{v}_{i}-\mathbf{v}_{j})-g(|\mathbf{v}_{i}|-v_{0})\dfrac{\mathbf{v}_{i}}{|\mathbf{v}_{i}|}+\bm{\xi}_{i}\
,$ (9)

where $\bm{\xi}_{i}$ is a Gaussian white noise with variance
$\langle\xi_{i}^{\alpha}(t)\
\xi_{j}^{\beta}(t^{\prime})\rangle=2T\delta_{ij}\delta_{\alpha\beta}\delta(t-t^{\prime})$.
We simulated a discretized version of these equations, using an Euler
discretization algorithm with $\Delta t=1.0\times 10^{-3}$, and computing at
each time step the adjacency matrix
$n_{ij}(t)=n_{ij}(|\mathbf{r}_{i}-\mathbf{r}_{j}|)$ according to a metric
rule,

$\displaystyle\begin{cases}n_{ij}=1\text{ }\text{if
}|\mathbf{r}_{i}-\mathbf{r}_{j}|\leq r_{c}\\\ n_{ij}=0\text{ }\text{if
}|\mathbf{r}_{i}-\mathbf{r}_{j}|>r_{c}\end{cases}$ (10)

All simulations have been performed with $r_{c}=1.0$.

###  A.1 Phase diagram

After a preliminary exploration of the space of parameters, we focused on a
few sets for which non-equilibrium effects are clearly visible in the whole
phase diagram. Most of the simulations discussed in the main text are
performed with $v_{0}=4$, $J=1.5$, $g=4$ and density $\rho=0.25$. In Fig. 3a,b
we report, for this set of parameters, the scalar polarization
$|{\bm{\Phi}}|=|(1/N)\sum_{i}\mathbf{v}_{i}/|\mathbf{v}_{i}|$ (i.e. the order
parameter) and its fluctuations as a function of temperature, for various
system sizes $N$. The fluctuations have been multiplied by a factor $N/T$ to
mimic what in equilibrium systems represents the static susceptibility. As it
can be seen from the figure, the model exhibits the typical behavior of the
Vicsek class, with a kinetic phase transition from a disordered state at large
noise to a state of collective motion at low temperatures. The peak of the
susceptibility at a given value of $N$ can be used to identify the critical
temperature $T_{c}(N)$ at that size.

To check the role of density, we also performed simulations at a larger value
of the density, with a second set of parameters: $v_{0}=4$, $J=0.5$, $g=4$ and
$\rho=1.5$ (the smaller value of $J$ ensures that the transition temperature
is not too large). The phase diagram for this set is displayed in Fig. 3c,d.
The qualitative picture remains the same, but with a larger transition
temperature.

Figure 3: Phase diagram of the model. a) Scalar polarization as a function of
temperature, for different values of the system size. b) Fluctuations of the
scalar polarization as a function of temperature, for different values of the
system size. The fluctuations have been multiplied by a factor $N/T$ to mimic
what in equilibrium systems represents the static susceptibility. The peak
indicates a kinetic transition at approximately $T_{c}=10$. The parameters of
the simulation are $v_{0}=4$, $J=1.5$, $g=4$ and $\rho=0.25$. Panels c) and
d): same as a) and b) for the second set of parameters: $v_{0}=4$, $J=0.5$,
$g=4$ and $\rho=1.5$. In this case the ordering transition occurs at a larger
temperature.

###  A.2 Density fluctuations and scaling

The behavior shown in Fig. 3 is consistent with a second-order transition
scenario. However, Vicsek models with metric interactions are known to
asymptotically exhibit a first-order transition, dominated by phase separation
Chaté  _et al._ (2008). Therefore, for large enough sizes we expect density
heterogeneities to become relevant and change the phenomenology at the
critical point. More precisely, it can be shown Di Carlo and Scandolo (2022)
that there is a crossover length-scale $\cal{R}$, whose value depends on the
microscopic parameters of the model, such that for $L<{\cal R}$ density
fluctuations are not strong enough and the large scale physics of the system
is described by the incompressible field theory Chen  _et al._ (2015), while
for $L>{\cal R}$ density fluctuations become relevant and the incompressible
regime does not hold anymore. Our numerical results indicate that - for the
sizes that we simulated - the system is still in the homogeneous -
incompressible regime.

Figure 4: a) Histogram of the number of interacting neighbors $n_{c}$, for
$\rho=0.25$ and $N=512$. b) Snapshot of a configuration in the stationary
regime. c) and d) same as a) and b) but with $\rho=0.5$ and $N=1728$.

To better investigate the role of heterogeneities we monitored the
distribution of the number of interacting neighbors $n_{i}$ during the
simulation. In Fig. 4a,b we report this quantity close to the critical
temperature for $\rho=0.25$ and $N=512$, the values used for most simulations,
together with a snapshot of a configuration in the stationary regime. The
histogram is well-peaked with small tails, indicating the absence of
heterogeneities, as also visually confirmed from panel b. For this value of
the density, heterogeneities start to be seen at the largest size considered
$N=1728$; however, they are not yet strong enough to alter the large scale
behavior. Indeed, as shown in the main text, dynamic scaling of the deviation
curves is very nicely obeyed with the critical exponent $z=1.7$ predicted by
the incompressible field theory. We also verified that the standard velocity
correlation functions obey dynamic scaling with the same exponent.

We note that $\rho=0.25$ corresponds to an average number of interacting
neighbors of just a few individuals (see Fig. 4a). We expect that increasing
the density, such number should also increase making density fluctuations
relatively less important. This is confirmed by simulations performed at
larger densities. For $\rho=1.5$ (corresponding to an average $n_{i}$ of $9$
individuals) the system remains quite homogeneous also at the largest size
considered, as reported in Fig. 4c,d. In Fig. 5 we show that - again - dynamic
scaling is very well obeyed with $z=1.7$. Our results indicate that increasing
the density makes the crossover length-scale ${\cal R}$ larger, meaning that a
second-order like scenario and the incompressible field theory fixed point
describe the large scale behavior up to larger system sizes.

Figure 5: a) Deviation curves $d(t)$ at $T=T_{c}(N)$ for several values of the
system size. b) Rescaled deviation $d(t)/d_{\rm max}(N)$ as a function of the
rescaled time $tL^{-z}$, with $z=1.7$. The parameters of the simulation are
$v_{0}=4,J=0.5,g=4,\rho=1.5$.

###  A.3 Computing correlations and responses above and below $T_{c}$

To compute correlations and response functions we followed the general
definitions given in the main text. As already mentioned, it is numerically
more convenient to look at the integrated response
$\chi_{\alpha\beta}(t)=\int_{0}^{t}d\tau R_{\alpha\beta}(\tau)$. In linear
response theory, $\chi$ is directly related to the change of the mean group
velocity under a constant perturbation, i.e. for
${\mathbf{h}}(t)={\mathbf{h}}\ \theta(t)$ one has
$V_{\alpha}(t)-V_{\alpha}(0)=\sum_{\beta}\chi_{\alpha\beta}(t)h_{\beta}$.
Therefore, it can be easily computed from data by monitoring the mean velocity
only. Since the probability distribution is symmetric under spacial
reflections, the off-diagonal components of the response (and correlation)
matrices are zero by symmetry. For $T>T_{c}$ the three diagonal components
$R_{\alpha\alpha}$ (and $\chi_{\alpha\alpha}$) are equal. However, for
$T<T_{c}$ the longitudinal component along the direction of the spontaneous
polarization vector is different from the perpendicular one (exactly as in
equilibrium alignment models), and they must be treated independently. In our
analysis we focus on the longitudinal response.

Measuring global correlations and responses in the ordered phase has
additional difficulties related to the finite size of the system and the so-
called wandering of the order parameter. Let us now discuss this issue. A key
feature of alignment models is that they exhibit, at low temperatures, a
spontaneous symmetry breaking in the thermodynamic limit. When the system
polarizes at the ordering transition, giving rise to a state of collective
motion, one specific direction is selected out of the many equivalent ones
consistent with the rotational symmetry of velocity alignment. A remnant of
this continuous symmetry for $T<T_{c}$ is the presence of soft modes (the so-
called Goldstone modes) in the manifold orthogonal to the global polarization,
which makes perpendicular fluctuations stronger than longitudinal ones Brézin
_et al._ (1973); Brézin and Wallace (1973). For finite sizes, ergodicity is
not strictly broken and the system can actually dynamically explore this
manifold, slowly changing the direction of the global polarization: this is
what is usually referred to as the ‘wandering’ of the order parameter. This
diffusive rotational motion occurs with diffusion coefficient of order $1/N$
and it becomes visible on numerical simulations on finite systems for large
enough times, making extremely difficult to compute reliably the bona fide
correlation functions of the vectorial order parameter.

To address this issue, for $T<T_{c}$ we explicitly break the symmetry with a
small external field $h_{0}$ to make sure the order parameter fluctuates
around this direction without wandering away. Then, we compute all the
relevant quantities (correlations, responses, deviation functions) in presence
of the field, and measure the off-equilibrium indicators. Finally, we verify
that the results are robust upon decreasing the field. We note that this
procedure (computing at finite $h$ and then sending the field to zero) is also
used in theoretical equilibrium studies of Heisenberg models Brézin  _et al._
(1973). More precisely, we performed simulations for two different values of
the external field $h_{0}=0.10$ and $h_{0}=0.05$. The external field
explicitly breaks the symmetry along the direction $\hat{x}$, and correlation
and response functions are computed for the $x$ component of the vectorial
order parameter. Fig. 6 shows the deviation from equilibrium encapsulated by
the effective temperature, as defined in the main text. The low-temperature
phase branch depends on the value of $h_{0}$. However, reducing $h_{0}$ lowers
the obtained effective temperature. Hence, the trend strongly suggests that,
in the (computationally hard) limit $h_{0}\to 0$ the peak at the transition
point becomes even sharper, thus supporting the discussion and the general
picture described in the main text.

Figure 6: Effective temperature $T_{\rm eff}(\infty)$ as a function of
temperature; the relative difference with respect to the equilibrium value is
reported on the $y$-axis. The parameters of the simulation are
$J=1.5,g=4,\rho=0.25,N=512$. The low temperature branch $T<T_{c}$ has been
computed with two different values of the external field.

Finally, we notice that the FDT violation seems to approach a non-zero limit
for $T\to 0$, contrary to what happens for the EPR (Fig 1)c). This behavior is
related to the fact that, while the EPR is computed in zero field, the FDT
violation below $T_{c}$ is evaluated (for the reasons discussed above) with
the small field $h_{0}$. The external field, in presence of periodic boundary
conditions, introduces an additional probability current and consequent FDT
violations, not explicitly considered in the EPR formula (see Fischer  _et
al._ (2018) for a similar case). However, this effect decreases with $h_{0}$
(as also shown in Fig. 6) and it therefore vanishes when $h_{0}\to 0$.

##  APPENDIX B Response function from the OM action

The Onsager-Machlup (OM) formalism is a well-known path-integral approach to
Langevin stochastic dynamics that allows computing expectation values starting
from the dynamical equations. In this section we briefly describe the method,
as applied to our case. More details can be found in standard textbooks of
stochastic processes.

Let us rewrite the equations of motion (Eq. (1) of the main text) in a more
compact form

$\dfrac{d\mathbf{v}_{i}}{dt}=\mathbf{F}_{i}(\mathbf{v},\mathbf{r})+\bm{\xi}_{i}\
\quad\quad\frac{d\mathbf{r}_{i}}{dt}=\mathbf{v}_{i}(t)$ (11)

where $\mathbf{F}_{i}=-\partial{\cal H}/\partial\mathbf{v}_{i}$,
$\mathbf{v}(t)\equiv\\{\mathbf{v}_{i}(t)\\}$,
$\mathbf{r}(t)\equiv\\{\mathbf{r}_{i}(t)\\}$, and
$\bm{\xi}(t)\equiv\\{\bm{\xi}_{i}(t)\\}$ \- we remind - is a Gaussian white
noise with variance
$\langle\xi_{i}^{\alpha}(t)\xi_{j}^{\beta}(t^{\prime})\rangle=2T\delta_{\alpha\beta}\delta_{ij}\delta(t-t^{\prime})$.
The main idea behind the OM approach is that any physical observable $O(t)$ at
time $t$ is a random variable, which can only depend on the noise realization
for times before $t$, and on the initial conditions. Its average value
$\langle O(t)\rangle$ is therefore obtained by averaging over the functional
distribution of the noise realizations. However, since observables are
ultimately functions of velocities and positions, such averages can be
conveniently re-expressed in terms of the distribution probability of the
velocities and positions themselves. Given a thermal history $\bm{\xi}(t)$,
velocities are obtained via Eq. 11 as
$\mathbf{v}_{i\bm{\xi}}(t)=\mathbf{v}_{i}(0)+\int_{0}^{t}dt^{\prime}[\bm{F}_{i}(\mathbf{v},\mathbf{r})+\bm{\xi}_{i}]$.
Positions are instead deterministically obtained from the velocities, i.e.
$\mathbf{r}_{i\mathbf{v}}(t)=\mathbf{r}_{i}(0)+\int_{0}^{t}dt^{\prime}\mathbf{v}_{i}(t^{\prime})$.
The distribution of velocities and positions, given the initial conditions,
can therefore be immediately derived as

$P\left[\mathbf{v}(t),\mathbf{r}(t)|\mathbf{v}(0),\mathbf{r}(0)\right]=\int
D\left[\bm{\xi}(t^{\prime})\right]P\left[\bm{\xi}(t^{\prime})\right]\prod_{i}\delta\left[\mathbf{v}_{i}(t^{\prime})-\mathbf{v}_{i\bm{\xi}}(t^{\prime})\right]\
\delta[\mathbf{r}_{i}(t)-\mathbf{r}_{i\mathbf{v}}(t)].$ (12)

where thermal histories are considered between times $0$ and $t$. Inverting
the dependence between $\mathbf{v}(t)$ and $\bm{\xi}(t)$, and exploiting the
explicit form of the Gaussian noise distribution, we get

$\displaystyle
P\left[\mathbf{v}(t),\mathbf{r}(t)|\mathbf{v}(0),\mathbf{r}(0)\right]=$
$\displaystyle
e^{-\large\int_{0}^{t}dt^{\prime}\left[\frac{1}{4T}\sum_{i}\left(\dot{\mathbf{v}}_{i}-\bm{F}_{i}(\mathbf{v},\mathbf{r})\right)^{2}+f[\mathbf{v},\mathbf{r}]\right]}\delta[\mathbf{r}_{i}(t)-\mathbf{r}_{i\mathbf{v}}(t)]$
(13) $\displaystyle=$ $\displaystyle
e^{-S[\mathbf{v}(t),\mathbf{r}(t)|\mathbf{v}(0),\mathbf{r}(0)]}\delta[\mathbf{r}_{i}(t)-\mathbf{r}_{i\mathbf{v}}(t)]\
.$

The term $S[\mathbf{v}(t)),\mathbf{r}(t)|\mathbf{v}(0),\mathbf{r}(0)]$ is
called “Onsager-Machlup” (OM) action. The function $f[\mathbf{v},\mathbf{r}]$
comes from the Jacobian factor $|{\rm det}(\delta\bm{\xi}/\delta\mathbf{v})|$
and its explicit expression depends on the stochastic integration scheme.
Using Ito’s rules of calculus (the more convenient ones for simulations) the
Jacobian is equal to $1$ and $f=0$. The expected value of any observable over
the thermal history can be computed as,

$\langle O(t)\rangle=\int
D\left[\bm{v}(t^{\prime})\right]D\left[\bm{r}(t^{\prime})\right]O(t^{\prime})e^{-S[\mathbf{v}(t^{\prime}),\mathbf{r}(t^{\prime})|\mathbf{v}(0),\mathbf{r}(0)]}\delta[\mathbf{r}_{i}(t^{\prime})-\mathbf{r}_{i\mathbf{v}}(t^{\prime})]\
,$ (14)

with

$S[\mathbf{v}(t),\mathbf{r}(t)|\mathbf{v}(0),\mathbf{r}(0)]=\frac{1}{4T}\sum_{i,\alpha}\int_{0}^{t}dt^{\prime}\
\left(\dot{v}_{i}^{\alpha}(t^{\prime})-F_{i}^{\alpha}(t^{\prime})\right)^{2}$
(15)

Expression (14) can be eventually averaged over the distribution of the
initial conditions. In the main text, for notational convenience, we omitted
the dependence on the initial conditions, and included the delta function on
the positions in the measure. We note however that in the stationary state, if
the system is ergodic, the dependence on initial conditions is lost.

###  B.1 The response function

Using Eq. (14) we can immediately compute the response function, i.e.

$\displaystyle R^{\alpha\beta}_{ij}(t,s)=$ $\displaystyle\dfrac{\delta\langle
v_{i}^{\alpha}(t)\rangle_{h}}{\delta
h_{j}^{\beta}(s)}\bigg{|}_{h=0}=-\Bigg{\langle}v_{i}^{\alpha}(t)\dfrac{\delta
S}{\delta h_{j}^{\beta}(s)}\Bigg{\rangle}\Bigg{|}_{h=0}$ (16) $\displaystyle=$
$\displaystyle\frac{1}{2T}\langle
v_{i}^{\alpha}(t)(\dot{v}_{j}^{\beta}(s)-F_{j}^{\beta}(s))\rangle=\frac{1}{2T}\dfrac{\partial
C_{ij}^{\alpha\beta}(t,s)}{\partial s}-\frac{1}{2T}\langle
v_{i}^{\alpha}(t)F_{j}^{\beta}(s))\rangle\ ,$

where, in the last line, averages are performed in zero field. In the
stationary state, the response and correlation function depend on $t$ and $s$
only through their difference. A simpler expression for the response function
can then be obtained. Let’s assume without loss of generality that $t>s$. In
this case, by causality, $R^{\alpha\beta}_{ij}(s-t)=0$. By subtracting it from
Eq. (16) we then get, for $t>s$,

$TR^{\alpha\beta}_{ij}(t-s)=-\dfrac{1}{2}\left[\dfrac{dC_{ij}^{\alpha\beta}(t-s)}{dt}+\dfrac{dC_{ij}^{\alpha\beta}(s-t)}{dt}\right]-d_{ij}^{\alpha\beta}(t-s)\
,$ (17)

where $d_{ij}^{\alpha\beta}(t-s)=(1/2)\langle
v_{i}^{\alpha}(t)F_{j}^{\beta}(s)-v_{i}^{\alpha}(s)F_{j}^{\beta}(t)\rangle$.
The off-diagonal terms with $\alpha\neq\beta$ are identically zero by symmetry
arguments. For the diagonal ones we get,

$TR^{\alpha\alpha}_{ij}(t-s)=-\dfrac{dC_{ij}^{\alpha\alpha}(t-s)}{dt}-d_{ij}^{\alpha\alpha}(t-s)\
,$ (18)

which is Eq. (3) of the main text. An explicit computation of the term
$d_{ij}^{\alpha\alpha}(t-s)$ for the model at hand can be obtained using the
expression of the force

${\mathbf{F}}_{i}=-J\sum_{j}n_{ij}({\mathbf{v}}_{i}-{\mathbf{v}}_{j})-g(|{\mathbf{v}}_{i}|-v_{0})\frac{{\mathbf{v}}_{i}}{|{\mathbf{v}}_{i}|}\
,$ (19)

and it gives,

$\displaystyle d_{ij}^{\alpha\alpha}(t-s)=$
$\displaystyle\dfrac{J}{2}\sum_{k}\langle
v_{i}^{\alpha}(s)n_{jk}(t)(v_{j}^{\alpha}(t)-v_{k}^{\alpha}(t))\rangle-\dfrac{J}{2}\sum_{k}\langle
v_{i}^{\alpha}(t)n_{jk}(s)(v_{j}^{\alpha}(s)-v_{k}^{\alpha}(s))\rangle$ (20)
$\displaystyle+\dfrac{gv_{0}}{2}\left(\langle
v_{i}^{\alpha}(t)\dfrac{v_{j}^{\alpha}(s)}{|\bm{v}_{j}(s)|}\rangle-\langle
v_{i}^{\alpha}(s)\dfrac{v_{j}^{\alpha}(t)}{|\bm{v}_{j}(t)|}\rangle\right)$

To obtain the global response to a homogeneous field (see main text), we sum
both sides of Eq. (18) over the $i$ and $j$ indexes and divide by $N$. We then
get,

$TR_{\alpha\alpha}(t-s)=-\frac{d}{dt}C_{\alpha\alpha}(t-s)-d_{\alpha\alpha}(t-s)$
(21)

where $C_{\alpha\alpha}(t-s)=N\langle V_{\alpha}(t)V_{\alpha}(s)\rangle$ is
the correlation of the global velocity and $d_{\alpha\alpha}(t-s)$ can be
computed directly from Eq. (20). Interestingly, given that the matrix $n_{ij}$
is symmetric, many terms simplify and we get

$d_{\alpha\alpha}(t-s)=\frac{1}{N}\sum_{ij}d_{ij}^{\alpha\alpha}(t-s)=\dfrac{gv_{0}N}{2}\left[\langle
V_{\alpha}(t)\Phi_{\alpha}(s)\rangle-\langle
V_{\alpha}(s)\Phi_{\alpha}(t)\rangle\right],$ (22)

recovering Eq. (6) of the main text.

##  APPENDIX C Entropy Production Rate

The average entropy production $\Sigma(t)$ can be defined as the averaged
logarithm of the ratio of the probability of observing a trajectory of length
$t$, and the probability of the reversed process Lebowitz and Spohn (1999);
Seifert (2012). For a Markov process, it can be conveniently expressed in
terms of the OM action,

$\displaystyle\Sigma(t)$
$\displaystyle=\Bigg{\langle}\log\dfrac{P[{\mathbf{v}}(t),\mathbf{r}(t)]}{P[{\mathbf{v}}^{\dagger}(t),\mathbf{r}^{\dagger}(t)]}\Bigg{\rangle}=\Bigg{\langle}\log\dfrac{P[{\mathbf{v}}(t),\mathbf{r}(t)|{\mathbf{v}}(0),\mathbf{r}(0)]P_{0}({\mathbf{v}}(0),\mathbf{r}(0))}{P[{\mathbf{v}}^{\dagger}(t),\mathbf{r}^{\dagger}(t)|{\mathbf{v}}^{\dagger}(0),\mathbf{r}^{\dagger}(0)]P_{0}({\mathbf{v}}^{\dagger}(0),\mathbf{r}^{\dagger}(0))}\Bigg{\rangle}$
(23) $\displaystyle=\langle
S[{\mathbf{v}}^{\dagger}(t),\mathbf{r}^{\dagger}(t)|{\mathbf{v}}^{\dagger}(0),\mathbf{r}^{\dagger}(0)]-S[{\mathbf{v}}(t),\mathbf{r}(t)|{\mathbf{v}}(0),\mathbf{r}(0)]\rangle+\Bigg{\langle}\log\dfrac{P_{0}({\mathbf{v}}(0),\mathbf{r}(0))}{P_{0}({\mathbf{v}}^{\dagger}(0),\mathbf{r}^{\dagger}(0))}\Bigg{\rangle}$

where averages are taken over the ensemble of the stochastic paths. The last
term in the r.h.s of Eq. (23) is a boundary term that does not depend on time:
since we are ultimately going to take the temporal average rate, we disregard
it in the following. The first two terms at the r.h.s. can be computed from
the explicit expression of the OM action. The state variables transform under
time reversal as
$\bm{r}(t^{\prime})\mapsto\bm{r}^{\dagger}(t^{\prime})=\bm{r}(t-t^{\prime})$
and
${\mathbf{v}}(t^{\prime})\mapsto{\mathbf{v}}^{\dagger}(t^{\prime})=-{\mathbf{v}}(t-t^{\prime})$.
For what concerns the forces ${\mathbf{F}_{i}(t^{\prime})}$ (see Eq.(19)),
given that positions enter only via the adjacency matrix $n_{ij}$ (which is
invariant under time reversal), we have
${\mathbf{F}_{i}(t^{\prime})}\mapsto{\mathbf{F}}_{i}^{\dagger}(t^{\prime})=-{\mathbf{F}_{i}(t-t^{\prime})}$.
From Eq.(15) we then get

$\displaystyle
S[{\mathbf{v}}^{\dagger}(t),\mathbf{r}^{\dagger}(t)|{\mathbf{v}}^{\dagger}(0),\mathbf{r}^{\dagger}(0)]-S[{\mathbf{v}}(t),\mathbf{r}(t)|{\mathbf{v}}(0),\mathbf{r}(0)]$
$\displaystyle=\frac{1}{4T}\sum_{i}\sum_{\alpha}\int_{0}^{t}dt\left[\left(\frac{dv_{i}^{\alpha}}{dt}+F_{i}^{\alpha}\right)^{2}-\left(\frac{dv_{i}^{\alpha}}{dt}-F_{i}^{\alpha}\right)^{2}\right]=$
(24) $\displaystyle=\frac{1}{T}\sum_{i}\sum_{\alpha}\int_{0}^{t}dt\
\frac{dv_{i}^{\alpha}}{dt}\circ F_{i}^{\alpha}\ ,$

where $\circ$ denotes the Stratonovich prescription, always arising regardless
of the type of stochastic calculus in the original integral (indeed the
inverse process always has by construction the complementary prescription of
the forward one).

We can now evaluate the average entropy production rate. In the stationary
state, and assuming ergodicity, we get

$\displaystyle{\rm EPR}$
$\displaystyle=\lim_{t\to\infty}\frac{\Sigma(t)}{t}=\lim_{t\to\infty}\frac{1}{t}\frac{1}{T}\sum_{i}\sum_{\alpha}\int_{0}^{t}dt\
\Big{\langle}\frac{dv_{i}^{\alpha}}{dt}\circ F_{i}^{\alpha}\Big{\rangle}$ (25)
$\displaystyle=\frac{1}{T}\sum_{i}\sum_{\alpha}\Big{\langle}\frac{dv_{i}^{\alpha}}{dt}\circ
F_{i}^{\alpha}\Big{\rangle}$
$\displaystyle=-\frac{1}{T}\sum_{ij}\sum_{\alpha}\Big{\langle}\frac{dv_{i}^{\alpha}}{dt}\circ\left[Jn_{ij}(v_{i}^{\alpha}-v_{j}^{\alpha})+g\frac{v_{i}^{\alpha}}{|\bm{v}_{i}|}(|\bm{v}_{i}|-v_{0})\right]\Big{\rangle}.$

This expression can be further manipulated observing that, in the stationary
state,

$0=\Big{\langle}\dfrac{d\mathcal{H}}{dt}\Big{\rangle}=\frac{J}{2}\sum_{i,j}\sum_{\alpha}\Big{\langle}\frac{dn_{ij}}{dt}\circ(v_{i}^{\alpha}-v_{j}^{\alpha})^{2}\Big{\rangle}-\sum_{i}\sum_{\alpha}\Big{\langle}F_{i}^{\alpha}\circ\dfrac{dv_{i}^{\alpha}}{dt}\Big{\rangle}.$
(26)

Hence, substituting in Eq. (25), we find

${\rm
EPR}=\frac{J}{2T}\sum_{i,j}\Big{\langle}\frac{dn_{ij}}{dt}\circ\sum_{\alpha}(v_{i}^{\alpha}-v_{j}^{\alpha})^{2}\Big{\rangle}.$
(27)

##  APPENDIX D Harada-Sasa relationships

###  D.1 Harada-Sasa relationship from the OM action

The Harada-Sasa relationship is a formula that connects the EPR with a
violation of the FDT, firstly derived for Langevin dynamics in Harada and Sasa
(2005). To illustrate such relation in our case, let us consider the
expression for the response function obtained in the previous section, and
take the difference $R^{\alpha\alpha}_{ii}(t-s)-R^{\alpha\alpha}_{ii}(t-s)$.
By causality only one of the two responses is different from zero and we get
(without assuming any order between the two times)

$R^{\alpha\alpha}_{ii}(t-s)-R^{\alpha\alpha}_{ii}(s-t)=-\dfrac{1}{T}\dfrac{dC_{ii}^{\alpha\alpha}(t-s)}{dt}-\dfrac{1}{T}d_{ii}^{\alpha\alpha}(t-s).$
(28)

We take now $s=t-\Delta t$, divide both sides by $\Delta t$, and take the
limit as $\Delta t$ goes to $0$. We consider first the left-hand side,

$\displaystyle\lim_{\Delta t\to 0}\dfrac{R^{\alpha\alpha}_{ii}(\Delta
t)-R^{\alpha\alpha}_{ii}(-\Delta t)}{\Delta t}=\lim_{\Delta t\to
0}\int_{-\infty}^{\infty}\dfrac{d\omega}{2\pi}\tilde{R}^{\alpha\alpha}_{ii}(\omega)\dfrac{e^{-i\omega\Delta
t}-e^{i\omega\Delta t}}{\Delta t}=$ (29)
$\displaystyle=-\int_{-\infty}^{\infty}\dfrac{d\omega}{2\pi}2i\omega\tilde{R}^{\alpha\alpha}_{ii}(\omega)=\int_{-\infty}^{\infty}\dfrac{d\omega}{2\pi}2\omega\
\mathfrak{Im}\left(\tilde{R}^{\alpha\alpha}_{ii}(\omega)\right),$

where
$\tilde{R}^{\alpha\alpha}_{ii}(\omega)=\int_{-\infty}^{\infty}dt\exp{(i\omega
t)}R^{\alpha\alpha}_{ii}(t)$ is the Fourier transform of the response, and
$\mathfrak{Im}$ denotes the imaginary part (the only one surviving inside the
integral, being the response a real function). The two terms at the r.h.s of
Eq. (28) can be handled in a similar way,

$\lim_{\Delta t\to 0}\frac{1}{\Delta
t}\dfrac{1}{T}\dfrac{dC_{ii}^{\alpha\alpha}(\Delta
t)}{dt}=-\dfrac{1}{T}\int_{-\infty}^{\infty}\dfrac{d\omega}{2\pi}\omega^{2}\tilde{C}_{ii}^{\alpha\alpha}(\omega),$
(30)

and

$\displaystyle\lim_{\Delta t\to 0}\frac{1}{\Delta
t}\dfrac{1}{T}d^{\alpha\alpha}_{ii}(\Delta t)=$
$\displaystyle\dfrac{1}{T}\lim_{\Delta t\to 0}\dfrac{\langle
v_{i}^{\alpha}(t)F_{i}^{\alpha}(t-\Delta t)\rangle-\langle
v_{i}^{\alpha}(t-\Delta t)F_{i}^{\alpha}(t)\rangle}{2\Delta t}=$ (31)
$\displaystyle=$ $\displaystyle\dfrac{1}{T}\lim_{\Delta t\to 0}\dfrac{\langle
v_{i}^{\alpha}(t)F_{i}^{\alpha}(t-\frac{\Delta t}{2})\rangle-\langle
v_{i}^{\alpha}(t-\frac{\Delta t}{2})F_{i}^{\alpha}(t)\rangle}{\Delta t}=$
$\displaystyle=$ $\displaystyle\dfrac{1}{T}\lim_{\Delta t\to 0}\dfrac{\langle
v_{i}^{\alpha}(t+\frac{\Delta t}{2})F_{i}^{\alpha}(t)\rangle-\langle
v_{i}^{\alpha}(t-\frac{\Delta t}{2})F_{i}^{\alpha}(t)\rangle}{\Delta t}=$
$\displaystyle=$ $\displaystyle\dfrac{1}{T}\langle\dot{v}_{i}^{\alpha}\circ
F_{i}^{\alpha}\rangle,$

where stationarity has been used, together with the explicit definition of the
Stratonovich differentiation. Putting all the pieces together, we find,

$\int_{-\infty}^{\infty}\dfrac{d\omega}{2\pi}2\omega\
\mathfrak{Im}\left(\tilde{R}^{\alpha\alpha}_{ii}(\omega)\right)=\dfrac{1}{T}\int_{-\infty}^{\infty}\dfrac{d\omega}{2\pi}\omega^{2}\tilde{C}_{ii}^{\alpha\alpha}(\omega)-\dfrac{1}{T}\langle\dot{v}_{i}^{\alpha}\circ
F_{i}^{\alpha}\rangle.$ (32)

Finally, comparing with Eq. (25), we get

${\rm
EPR}=\sum_{i,\alpha}\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\frac{\omega}{T}\left[\omega\tilde{C}_{ii}^{\alpha\alpha}(\omega)-2T\mathfrak{Im}(\tilde{R}_{ii}^{\alpha\alpha}(\omega))\right]\
,$ (33)

which is the Harada-Sasa formula, Eq. (8) of the main text.

###  D.2 Numerical check of the Harada-Sasa formula

In the previous sections we derived two distinct expressions of the EPR. The
first one, detailed in Sec.C, can be used to compute the EPR directly from the
trajectories from the numerical simulation, by implementing a discrete version
of the time integral appearing in Eq. (25). This method is very efficient, it
has been already tested in Ferretti  _et al._ (2022), and it produces the red
curve in Fig. 2b of the main text. The Harada-Sasa relation, on the other
hand, is more complicated to check numerically, especially for strongly
interacting systems as the one considered in this work. In this section we
discuss why, and how we addressed the various technical issues.

To compute numerically the r.h.s. of Eq. (33), one needs first of all to
compute the local responses $R^{\alpha\alpha}_{ii}(t)$. As discussed in the
main text, local responses are more demanding to calculate than global ones.
Luckily, one can resort to the so-called Malliavin weight sampling Warren and
Allen (2013). The basic idea of this method can be grasped by looking back at
the derivation of Eq. (16). From the second line, we see that the response can
be expressed as $R^{\alpha\alpha}_{ii}(t-s)=1/(2T)\langle
v_{i}^{\alpha}(t)\xi_{i}^{\alpha}(s)\rangle$, where averages are performed in
zero field. Therefore, by keeping track of the noise stories, one can compute
local responses without perturbing the system. We refer the reader to the
relevant literature for smart implementations of the method.

One issue, specifically related to polar systems of finite size, is the
spurious effect due to the wandering of the order parameter below the ordering
transition. To contain the effect of wandering, we notice that FDT relations
and the Harada-Sasa formula only involves derivatives of the correlation
function, and they therefore remain the same either considering simple or
connected correlations. At low temperatures we then used connected correlation
functions,

$C_{ii}^{\alpha\alpha}(t-s)=\langle\delta v_{i}^{\alpha}(t)\delta
v_{i}^{\alpha}(s)\rangle,\quad\quad\delta
v_{i}^{\alpha}(t)=v_{i}^{\alpha}(t)-\dfrac{1}{N}\sum_{j}v_{j}^{\alpha}(t)\ ,$
(34)

where the space average is used as a proxy of the average mean velocity. In
the thermodynamic limit, this procedure exactly recovers the theoretical
prediction Eq. (33). At the same time, by automatically subtracting the
wandering, it also provides a better estimate of this expression also for
finite systems. For consistency, and to ensure a better comparison, we also
evaluate the EPR from Eq. (25) using velocities in the center of mass
reference frame.

The Harada-Sasa relation involves the Fourier transforms of the correlation
and response functions, and the multiplication of their difference by a
(diverging) factor $\omega$. This implies that, whenever there is any small
error in determining precisely these function, the integral of Eq. (33) can
differ significantly from the theoretical value, possibly even becoming
divergent. The first trivial source of error is the finite time window of
observation of the correlation and response functions, affecting the discrete
resolution of the $\omega$ points, and thus the accuracy of the integral. Less
trivial, instead, is the impact of the finite integration step $\Delta t$ of
the numerical simulation: it is not commonly known that, even at equilibrium,
the temperature defined as the noise variance coefficient and temperature
defined as the factor appearing in the Gibbs-Boltzmann distribution are not
the same for discrete-time dynamics, but they differ by a term of $O(\Delta
t)$ Cavagna  _et al._ (2014). As a consequence, even at equilibrium, discrete
dynamics response differs from the continuum expression by a term $O(\Delta
t)$, which can affect even quite significantly the result of the integral. To
address this problem, we hence performed simulations with a smaller $\Delta
t=1.0\times 10^{-4}$, at the price of longer times needed for thermalization.
Finally, the finite time resolution of the discrete points of the correlation
functions imposes, in Fourier’s space, a constraint on the tails of the
transformed functions, making them hit zero at a certain cutoff; to avoid this
problem, a high resolution in time is needed for the correlation and response
functions, at the price of stronger noise when taking the derivative of the
correlation function and the integrated response, because of smaller steps.
All things considered, the signal is very noisy and it requires a lot of
samples. The way we estimated the integral at the r.h.s. of Eq. (33), and
produced the blue curve of Fig. 2b, is the following: we performed many runs
of numerical simulations, divided them into groups and computed the average
correlation and response functions for each group. Then, for each group
average, we computed a cut-offed version of the integral,

$F(\Omega)=\sum_{i,\alpha}\int_{-\Omega}^{\Omega}\frac{d\omega}{2\pi}\frac{\omega}{T}\left[\omega\tilde{C}_{ii}^{\alpha\alpha}(\omega)-2T\mathfrak{Im}(\tilde{R}_{ii}^{\alpha\alpha}(\omega))\
.\right]$ (35)

This function should have, was not for the noise, an asymptote for
$\Omega\to\infty$. We then assign a value to the integral by performing a
numerical fit of the asymptote. We repeat the procedure for each group and we
then take the average of them all. The error bars are computed as the standard
deviation across groups and, considered the many sources of error, they are
likely to be still underestimated. We also tried a different kind of analysis,
where to reduce the impact of noise we first fit the integrand with
appropriate functional forms, and then compute the integral. The results are
in general not very different for the average, while the associated error bars
are smaller, while realistically they should be bigger. Finally, by comparing
with the EPR computed from trajectories (see Fig. 2b), we see that there is
still a slight, but systematic, bias below $T_{c}$ where the EPR via the
Harada-Sasa formula is underestimated. We were not able to pinpoint the reason
of this. Nevertheless, in light of all the aforementioned numerical
challenges, we consider our results a quite satisfactory validation of the
formula.

##  APPENDIX E Equilibrium limits

###  E.1 Recovering equilibrium in the active frame

In the main text we derived an expression, Eq. (6), where FDT violations are
quantified in terms of simple correlation functions between the mean velocity
and the polarization. One way to read Eq. (6) is that, out of equilibrium, the
information content of the polarization is larger than the one of the
collective velocity, the first influencing the future behavior of the second
more than viceversa. In general, we expect violation of detailed balance to
manifest differently in different degrees of freedom: by introducing
fluctuating speeds we exploited this property to find easily accessible
signatures of FDT violations.

Interestingly, in our system the slow critical modes (i.e. the polarization)
are the ones most affected by detailed balance violations. On the contrary,
the speeds - which are fast modes and mark the difference between
${\mathbf{\Phi}}$ and ${\bf V}$, are not. To illustrate this, we now consider
the behavior of the mean velocity in a modified frame, i.e. when subtracting
the mean “active trajectory” $v_{0}{\mathbf{\Phi}}$ from the individual
velocities. We therefore define the reduced mean velocity ${\tilde{\bf
V}}={\bf V}-v_{0}{\mathbf{\Phi}}$ and show that it responds as an equilibrium
quantity when we perturb the system with an external field conjugated to the
velocity. To do so, we first observe that, from the equation of motion of the
single particles, an equation of motion for the global velocity can be
derived,

$\frac{dV_{\alpha}(t)}{dt}=-gV_{\alpha}(t)+gv_{0}\Phi_{\alpha}(t)+\zeta_{\alpha}(t),$
(36)

with
$\langle\zeta_{\alpha}(t)\zeta_{\beta}(t^{\prime})\rangle=(2T/N)\delta_{\alpha\beta}\delta(t-t^{\prime})$.
Starting from this equation we can compute the antisymmetric average
$\langle\Phi_{\alpha}(s)dV_{\alpha}(t)/dt\rangle-\langle\Phi_{\alpha}(t)dV_{\alpha}(s)/ds\rangle$.
In the stationary state, we get (for $t>s$)

$2\Big{\langle}\Phi_{\alpha}(s)\frac{dV_{\alpha}(t)}{dt}\Big{\rangle}=-g\big{[}\langle\Phi_{\alpha}(s)V_{\alpha}(t)\rangle-\langle\Phi_{\alpha}(t)V_{\alpha}(s)\rangle\big{]}+\langle\Phi_{\alpha}(t)\zeta_{\alpha}(s)\rangle\
.$ (37)

The first term in the r.h.s. is proportional to the deviation function and it
can be re-expressed in terms of the correlation and the response of the mean
velocity (see Eqs. (6) and (18)). The second term is nothing else than the
response of the polarization, i.e.
$2T\delta\langle\Phi_{\alpha}(t)\rangle/\delta h_{\alpha}(s)|_{h=0}$, as it
can be deduced from Eq. (16). Reordering the terms, we therefore get

$\frac{\delta\langle\tilde{V}_{\alpha}(t)\rangle}{\delta
h_{\alpha}(s)}\Bigg{|}_{h=0}=\dfrac{N}{T}\langle\tilde{V}_{\alpha}(t)\dfrac{d}{ds}V_{\alpha}(s)\rangle\
,$ (38)

which is precisely the relation holding for equilibrium systems. We note that
this result reminds the restoring of equilibrium relations in the Lagrangian
frame described in systems of non-equilibrium colloidal and granular particles
Speck and Seifert (2006); Chetrite  _et al._ (2008); Crisanti  _et al._ (2012)

###  E.2 Pointers limit

One can ask whether an equilibrium limit exists for the model of Eqs.(1)(2).
Since motility is the primary source of violation of detailed balance, one
could think that merely sending $v_{0}$ to zero should recover equilibrium,
and that $d(t)$ trivially vanishes because of the $v_{0}$ prefactor (see Eq.
(6)). However, a finite speed not only determines the rearrangement of the
interaction network (and the consequent entropy production), but also the
length of the individual directional degrees of freedom. If we just send
$v_{0}\to 0$, particles slow down, but all velocities become smaller while
their confining potential tends to a Gaussian theory, which would become
problematic in terms of ordering. The limiting model would have a very
different phenomenology than the original one. The proper way to define an
equilibrium limit is more subtle, and requires focusing on rescaled variables.
Let us define

$\bm{\sigma}_{i}=\dfrac{\mathbf{v}_{i}}{v_{0}}.$ (39)

In these new variables, the equations of motion become,

$\frac{d\bm{\sigma}_{i}}{dt}=J\sum_{j}n_{ij}(\bm{\sigma}_{i}-\bm{\sigma}_{j})-g(|\bm{\sigma}_{i}|-1)\dfrac{\bm{\sigma}_{i}}{|\bm{\sigma}_{i}|}+\bm{\tilde{\xi}}_{i}\
,\quad\quad\frac{d\mathbf{r}_{i}}{dt}=v_{0}\bm{\sigma}_{i}$ (40)

with
$\langle\bm{\tilde{\xi}}_{i}(t)\cdot\bm{\tilde{\xi}}_{j}(t^{\prime})\rangle=2d(T/v_{0}^{2})\delta_{ij}\delta(t-t^{\prime})$.
If we now rescale the temperature, i.e. $T\to{\tilde{T}}=T/v_{0}^{2}$ (which
amounts to considering a random force in the original equation that, as the
deterministic ones, scales with the speed), the only remaining dependence on
$v_{0}$ is in the equation of motion for the positions and, implicitly, in the
interaction network $n_{ij}$. When taking the limit $v_{0}\to 0$, the only
thing happening is that positions ‘freeze’ and the network becomes fixed. The
equation for the $\\{\bm{\sigma}_{i}\\}$ then becomes the relaxational
dynamics of a Heisenberg model with a soft constraint on the spins, thus
recovering an equilibrium model. Looking at the off-equilibrium deviation
$d(t)$ written in terms of the rescaled variables, it reads

$d(t-s)=d_{\alpha\alpha}(t-s)=\frac{gN}{\tilde{T}}[\langle\sigma_{\alpha}(t)\Phi_{\alpha}(s)\rangle-\langle\sigma_{\alpha}(s)\Phi_{\alpha}(t)\rangle]=0\
$

where $\sigma_{\alpha}=(1/N)\sum_{i}\sigma_{i}^{\alpha}$. We note that now
there is no $v_{0}$ prefactor, nevertheless when $v_{0}\to 0$ the deviation
vanishes because the two correlation functions in the r.h.s. become identical
due to the restoring of detailed balance. When considering the EPR, expressed
in terms of the rescaled variables, when $v_{0}\to 0$ it vanishes trivially
because the time derivative of the adjacency matrix becomes zero (see Eq.
(27)).


11institutetext: AGH University of Science and Technology, Faculty of Physics
and Applied Computer Science, aleja Mickiewicza 30, 30-059 Cracow, Poland
22institutetext: Université Paris Saclay, CNRS, CEA, Institut de physique
théorique, 91191 Gif-sur-Yvette, France 33institutetext: Department of
Chemistry, Stony Brook University, Stony Brook, NY 11794, USA 44institutetext:
Physics Department, Brookhaven National Laboratory, Upton, NY 11976, USA
55institutetext: Instituto de Física, Universidade de São Paulo, Rua do Matão,
1371, Butantã, 05508-090, São Paulo, Brazil

# Transverse momentum fluctuation in ultra-central Pb+Pb collision at the LHC

Rupam Samanta 1122<EMAIL_ADDRESS>Somadutta Bhatta 33 Jiangyong Jia 3344
Matthew Luzum 55 Jean-Yves Ollitrault 22

###### Abstract

The ATLAS collaboration has recently observed that the variance of the
transverse momentum per particle ($[p_{t}]$), when measured as a function of
the collision multiplicity ($N_{ch}$) in Pb+Pb collisions, decreases by a
factor $2$ for the largest values of $N_{ch}$, corresponding to ultra-central
collisions. We show that this phenomenon is naturally explained by invoking
impact parameter ($b$) fluctuations, which contribute to the variance, and
gradually disappear in ultra-central collisions. It implies that $N_{ch}$ and
$[p_{t}]$ are strongly correlated at fixed $b$, which is explained by the
local thermalization of the QGP medium.

##  1 Introduction

The left panel of Fig. 1 displays the variance of transverse momentum per
particle ($[p_{t}]$) as a function of the collision multiplicity in Pb+Pb
collisions at the LHC ATLAS:2022dov . (It is the dynamical contribution to the
variance, after subtracting the trivial statistical fluctuation from the
finite multiplicity.) It decreases over a narrow range of multiplicity around
$N_{ch}=3700$. This phenomenon is not explained by models such as HIJING,
which predicts a smooth $1/N_{ch}$ dependence of the variance for all ranges
of $N_{ch}$ Bhatta:2021qfk . We show that it is naturally explained if the
medium produced during the collison thermalizes.

Figure 1: Left: ATLAS data (light blue symbols) and our model fit presented in
Sec.3 (solid red curve). The green curve denotes the contribution due to $b$
fluctuation and the blue dashed line denotes the contribution of intrinsic
fluctuation. The vertical dashed line is the knee of the multiplicity
distribution. Right : Two dimensional distribution of $[p_{t}]$ and $N_{ch}$
from simulation of Pb+Pb collisions at 5.02 TeV and b=0. Blue symbols are
hydrodynamic calculations. Instead of $[p_{t}]$ and $N_{ch}$, we plot the
differences $N_{ch}-\overline{N_{ch}}$ and $\delta
p_{t}=[p_{t}]-\overline{p_{t}}$ where $\overline{N_{ch}}=6692$ and
$\overline{p_{t}}=1074MeV/c$ are the average values. The blue curve is the 99%
confidence ellipse of a Gaussian distribution fitted to these symbols. The
orange curve is a similar 99% confidence ellipse for 1.4 $\times$ 106
collisions simulated with HIJING. The straight lines represent the average
values $\delta p_{t}(N_{ch},b=0)$.

##  2 Hydro vs HIJING results

The impact parameter of the collision, $b$, is not measured, and varies event
by event. Fluctuations observed in a sample of collisions event originate from
fluctuations of $b$, and from fluctuations at fixed $b$. In order to
understand their respective roles, we run simulations of Pb+Pb collisions at
5.02 TeV at fixed $b$, more specifically at $b=0$.

First, we simulate 150 events using relativistic viscous hydrodynamics
Bozek:2021mov  and compute $N_{ch}$ and $[p_{t}]$ in each event. Fig. 1
displays the distribution of the two quantities from hydrodynamic simulation
where thermalization of the QGP is assumed. It is clearly visible that both
$[p_{t}]$ and $N_{ch}$ have significant fluctuation. As the impact parameter
is fixed, the origin of these fluctuations is purely quantum, due to genuine
intrinsic fluctuation at the initial state of the collision.

Most importantly, there is a strong positive correlation between $[p_{t}]$ and
$N_{ch}$, whose physical origin is clear: Larger $N_{ch}$ means larger density
$N_{ch}/V$, as the collision volume, which is mostly defined by the impact
parameter of the collision, is fixed. A larger density results in a larger
initial temperature if the system thermalizes. This in turn implies a higher
energy per particle in the final state, and eventually a larger transverse
momentum per particle $[p_{t}]$, generating the positive correlation between
$[p_{t}]$ and $N_{ch}$.

Second, we run simulations using the HIJING model Gyulassy:1994ew , which
considers the collision as a superposition of independent nucleon-nucleon
collision, and in which the particles have no interaction between each other
after production. The striking difference in comparison to the hydro results
is that the correlation between $N_{ch}$ and $[p_{t}]$ is smaller by a factor
$\sim$ 10\. In HIJING, a larger density does not imply a significant increase
of $[p_{t}]$, which can be ascribed to the lack of thermalization.

##  3 Model fit

Figure 2: Left : Charged particle multiplicity ($N_{ch}$) distribution in
Pb+Pb collision at 5.02 TeV measured by the ATLAS collaboration, denoted by
the light blue symbols. The solid red line represents the model fit. The thin
blue lines corresponds to the Gaussian distributions $P(N_{ch}|b)$ drawn for 0
%, 5 %, 10 % and 15 % centrality fractions. The vertical dashed line
corresponds to the "knee" ($\overline{N_{ch}{b=0}}$) of the distribution.
Right: Two dimensional joint distribution of transverse momentum per particle
$[p_{t}]$ and multiplicity $N_{ch}$ returned by the model fit to the
fluctuation data. Rather than plotting $[p_{t}]$, we plot $\delta
p_{t}=[p_{t}]-\overline{p_{t}}(0)$. The white curves represent the 99 %
confidence ellipses at fixed impact parameters. Corresponding to these values
$b$, a schematic representation of the colliding nuclei is shown. The solid
black line denotes the mean value of $\delta p_{t}$ averaged over events and
the red band represents 1-$\sigma$ band.

We model fluctuations of $[p_{t}]$ by assuming that the joint probability
distribution of $[p_{t}]$ and $N_{ch}$ at fixed $b$ is a two-dimensional
correlated Gaussian Samanta:2023amp , denoted by $P([p_{t}],N_{ch}|b)$. It is
characterized by five parameters: mean and variance of $[p_{t}]$ and $N_{ch}$,
denoted by $\overline{p_{t}}(b)$, $\overline{N_{ch}}(b)$, Var($p_{t}|b$),
Var($N_{ch}|b$) and the Pearson correlation coefficient between $[p_{t}]$ and
$N_{ch}$, $r_{N_{ch}}(b)$ which we expect to be positive in hydrodynamics, as
shown in Fig. 1 (right).

The mean and variance of $N_{ch}$ are obtained by fitting the measured
distribution of $N_{ch}$ as a superposition of Gaussians, along the lines of
Ref. Das:2017ned . The average value of $N_{ch}$ at $b=0$ is the knee of the
distribution, which is reconstructed precisely. Note that the fall of the
variance in Fig. 1 precisely occurs around the knee.

The mean transverse momentum $\overline{p_{t}}(b)$ is essentially independent
of $b$ for 30 % most central collisions ALICE:2018hza . We subtract this
constant and only study the distribution of the distribution of $\delta p_{t}$
= $[p_{t}]-\overline{p_{t}}$. We assume that the variance Var($p_{t}|b$)
varies with the mean according to a power law, Var($p_{t}|b$) =
$\sigma^{2}_{\delta
p_{t}}(\overline{N_{ch}}(0)/\overline{N_{ch}}(b))^{\alpha}$, where $\alpha$
and $\sigma_{\delta p_{t}}$ are constants Samanta:2023amp . For simplicity, we
assume the correlation $r_{N_{ch}}$ to be independent of $b$.

The variance of $\delta p_{t}$ at fixed $N_{ch}$ is then evaluated in the
following way. We first obtain the distribution of $\delta p_{t}$ at fixed
$N_{ch}$ and $b$ using $P(\delta p_{t}|N_{ch},b)=P(\delta
p_{t},N_{ch}|b)/P(N_{ch}|b)$. It is again a Gaussian distribution, therefore,
it is fully characterized by its mean and variance, which we denote by
$\kappa_{1}$ and $\kappa_{2}$. Both are functions of $N_{ch}$ and $b$. The
crucial point is that the mean $\kappa_{1}$ is proportional to the
multiplicity fluctuation $N_{ch}-\overline{N_{ch}(b)}$. It is proportional to
the Pearson correlation coefficient $r_{N_{ch}}$ between $\delta p_{t}$ and
$N_{ch}$ (straight blue line in Fig. 1 right).

The last step is to average over $b$. For this, one needs the probability
distribution of $b$ at fixed $N_{ch}$, which is obtained using Bayes’ theorem
$P(b|N_{ch})p(N_{ch})=P(N_{ch}|b)P(b)$ along the lines of Ref. Das:2017ned . A
simple calculation then shows that

$\displaystyle\text{mean}\ \langle\delta p_{t}\rangle$
$\displaystyle=\langle\kappa_{1}\rangle_{b}$ (1)
$\displaystyle\text{Var}(p_{t})$
$\displaystyle=(\langle\kappa_{1}^{2}\rangle_{b}-\langle\kappa_{1}\rangle_{b}^{2})+\langle\kappa_{2}\rangle_{b}$

where angular brackets $\langle\dots\rangle_{b}$ denote the average over $b$
for fixed $N_{ch}$. Two terms contribute to the variance: The first term
originates from fluctuations of the impact parameter $b$ at fixed $N_{ch}$,
and the second term is the true intrinsic (quantum) fluctuation, which is not
the by-product of $b$ fluctuation.

The expression in Eq. (1) is finally fitted to ATLAS data using the three
parameters : $\sigma_{\delta p_{t}}$, $r_{N_{ch}}$ and $\alpha$. Fig. 1
displays the model fit to the data (red lines) and the two separate
contributions to Var($p_{t}|N_{ch}$). Our model precisely explains the steep
decrease of the variance around the knee. Below the knee, both contributions
have comparable magnitudes, but above the knee, the first term, due to impact
parameter fluctuations, gradually disappears, causing the steep fall in the
variance data. This occurs because the effect of $b$ fluctuation becomes
negligible in ultra-central region due to the strict lower limit of $b$ ($\geq
0$), where the distribution of $\delta p_{t}$ becomes a truncated Gaussian
Das:2017ned ; Samanta:2023kfk .

Fig. 2 displays the two-dimensional distribution of $[p_{t}]$ and $N_{ch}$
returned by our fit. The white curves represent 99% confidence ellipses at
fixed impact parameters and they are tilted with respect to the horizontal
axis as in the hydrodynamic calculation in Fig. 1. The tilts denote the
positive correlation between $[p_{t}]$ and $N_{ch}$, characterized by
$r_{N_{ch}}$. From this plot, it can be seen that the width of $[p_{t}]$
distribution is partly because several ellipses contribute at a given $N_{ch}$
(first term in Eq. 1) and the other part comes from the vertical width of a
single ellipse (second term in Eq. 1). As a corollary, we also predict a
slight increase of the mean transverse momentum $\langle\delta p_{t}\rangle$
with $N_{ch}$, denoted by the black line in Fig. 2.

##  4 Physical significance and thermodynamic interpretation

Our results imply that there is a strong correlation between $[p_{t}]$ and
$N_{ch}$ at fixed impact parameter, which is quantified by the Pearson
correlation coefficient $r_{N_{ch}}=0.676$. As explained in Sec. 2, such a
correlation is naturally present in a hydrodynamic description. Hence, the
recent ATLAS data further support the common hypothesis that local
thermalization is achieved in Pb+Pb collisions. It is the first such evidence
which does not involve anisotropic flow.

Additionally, our analysis and methodology pave the way to separate the
classical (geometrical) and quantum (intrinsic) fluctuation. Moreover, the
slight increase of mean transverse momentum with multiplicity could be used to
extract the speed of sound in QGP medium Gardim:2019brr , which has recently
been measured by the CMS collaboration with great precision CMS-PAS-HIN-23-003
.

## Acknowledgments

R. S. is supported by the Polish National Science Center under grant NAWA
PRELUDIUM BIS: PPN/STA/2021/1/00040/U/00001 and PRELUDIUM BIS:
2019/35/O/ST2/00357. S.B and J.J are supported by DOE DE-SC0024602. M. L.
thanks the São Paulo Research Foundation (FAPESP) for support under grants
2021/08465-9, 2018/24720-6, and 2017/05685-2, as well as the support of the
Brazilian National Council for Scientific and Technological Development
(CNPq).

## References

  * (1) G. Aad et al. (ATLAS), Phys. Rev. C 107, 054910 (2023), 2205.00039
  * (2) S. Bhatta, C. Zhang, J. Jia, Phys. Rev. C 105, 024904 (2022), 2112.03397
  * (3) P. Bozek, R. Samanta, Phys. Rev. C 105, 034904 (2022), 2109.07781
  * (4) M. Gyulassy, X.N. Wang, Comput. Phys. Commun. 83, 307 (1994), nucl-th/9502021
  * (5) R. Samanta, S. Bhatta, J. Jia, M. Luzum, J.Y. Ollitrault (2023), 2303.15323
  * (6) S.J. Das, G. Giacalone, P.A. Monard, J.Y. Ollitrault, Phys. Rev. C 97, 014905 (2018), 1708.00081
  * (7) S. Acharya et al. (ALICE), Phys. Lett. B 788, 166 (2019), 1805.04399
  * (8) R. Samanta, J.a.P. Picchetti, M. Luzum, J.Y. Ollitrault, Phys. Rev. C 108, 024908 (2023), 2306.09294
  * (9) F.G. Gardim, G. Giacalone, J.Y. Ollitrault, Phys. Lett. B 809, 135749 (2020), 1909.11609
  * (10) Tech. rep., CMS preliminary, presented in Quark Matter (2023), https://cds.cern.ch/record/2870141