# Intermittent relaxation and avalanches in extremely persistent active matter

Yann-Edwin Keta  Laboratoire Charles Coulomb (L2C), Université de Montpellier,
CNRS, 34095 Montpellier, France  Rituparno Mandal  Institute for Theoretical
Physics, Georg-August-Universität Göttingen, 37077 Göttingen, Germany  Peter
Sollich  Institute for Theoretical Physics, Georg-August-Universität
Göttingen, 37077 Göttingen, Germany  Department of Mathematics, King’s College
London, London WC2R 2LS, UK  Robert L. Jack  Yusuf Hamied Department of
Chemistry, University of Cambridge, Lensfield Road, Cambridge CB2 1EW, United
Kingdom  Department of Applied Mathematics and Theoretical Physics, University
of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom  Ludovic
Berthier  Laboratoire Charles Coulomb (L2C), Université de Montpellier, CNRS,
34095 Montpellier, France  Yusuf Hamied Department of Chemistry, University of
Cambridge, Lensfield Road, Cambridge CB2 1EW, United Kingdom

###### Abstract

We use numerical simulations to study the dynamics of dense assemblies of
self-propelled particles in the limit of extremely large, but finite,
persistence times. In this limit, the system evolves intermittently between
mechanical equilibria where active forces balance interparticle interactions.
We develop an efficient numerical strategy allowing us to resolve the
statistical properties of elastic and plastic relaxation events caused by
activity-driven fluctuations. The system relaxes via a succession of scale-
free elastic events and broadly distributed plastic events that both depend on
the system size. Correlations between plastic events lead to emergent dynamic
facilitation and heterogeneous relaxation dynamics. Our results show that
dynamical behaviour in extremely persistent active systems is qualitatively
similar to that of sheared amorphous solids, yet with some important
differences.

##  I Introduction

Active matter systems continue to surprise and challenge us with their range
of diverse behaviour, generating new insights into systems that are relevant
for both biology and fundamental statistical mechanics vicsek95 ; marchetti13
; volpe16 ; solon15 . In dense disordered systems with slow dynamics, systems
of active particles exhibit glassy behaviour reminiscent of passive
equilibrium systems, as well as interesting new phenomena berthierkurchan13 ;
berthier14 ; ni13 ; szamel2015glassy ; mandal16 ; leomach19 ;
berthier2019glassy ; janssen2019active ; mandal2020extreme ; ManSol21 ;
keta2022disordered .

A fascinating idea is that collective motion observed on large length scales
may be triggered by driving forces acting at the particle scale. This
represents a form of dynamic self-organisation that emerges out of equilibrium
from many-body effects. A wide variety of collective behaviour is observed in
active systems, from bacterial turbulence dunkel2013fluid  to flocking motion
vicsek2012collective . These effects are revealed by the emergence of non-
trivial correlations cavagna2018physics , either in the instantaneous velocity
field or in displacement fields over larger time lags.

In systems of self-propelled particles without aligning interactions,
sometimes known as scalar active matter, the difference between active and
passive systems appears through the persistence time $\tau_{p}$ of the self-
propulsion berthier14 ; szamel2014self ; fodor2016how . The passive behaviour
can be recovered in the non-persistent limit $\tau_{p}\to 0$. The other
extreme of very large values of $\tau_{p}$ leads to distinct behaviour, such
as jamming and intermittent relaxation mandal2020extreme ; mandal2021how ;
keta2022disordered , reminiscent of sheared athermal systems
maloney2006amorphous . In such persistent systems, collective motion emerges
when crowding effects at large density compete with persistent activity
szamel2015glassy ; flenner2016nonequilibrium ; berthier2017active ;
henkes2020dense ; keta2022disordered . This happens in both fluid and glassy
states, and increasing the persistence time typically reinforces these effects
mandal2020extreme ; henkes2020dense ; szamel2021long  because the correlation
length of velocity correlations grows with $\tau_{p}$. However, computational
studies of this regime are arduous because of the separation of timescales
between $\tau_{p}$ and emergent correlation times on the one hand, and
molecular timescales on the other.

For very large but finite $\tau_{p}$, it has recently been shown that this
difficulty can be tackled using a strategy known as activity-driven dynamics
(ADD) mandal2021how ; mandal2020multiple . For systems without thermal noise,
this provides direct access to the limit of very large $\tau_{p}$, and allows
simulations of all dynamical relaxation processes taking place on timescales
comparable with and greater than $\tau_{p}$. The central idea is that for
times $t\ll\tau_{p}$, the system evolves as if the propulsion forces were
fixed: if these forces are not too strong, then the dynamics will converge to
a force-balanced configuration in a time $\tau_{0}\ll\tau_{p}$. For larger
forces, the system yields and there is no convergence to any mechanical
equilibrium liao2018criticality ; mandal2020extreme ; villarroel2021critical ;
ADD is restricted to forces below that yield point. After reaching a force-
balanced state, the system configuration in ADD remains almost constant until
the propulsion forces change significantly, which requires a time of order
$\tau_{p}$. Eventually, these forces will change sufficiently to destabilise
the force-balanced state, at which point the system must switch to a new
force-balanced configuration, which involves significant particle motion.

The resulting intermittent dynamical motion has several important
consequences. First, it means that standard simulation methods become
inefficient because they require a time step $\delta
t\lesssim\tau_{0}\ll\tau_{p}$. ADD circumvents this problem by replacing the
explicit dynamical integration over times of order $\tau_{0}$ with an energy
minimisation step. This offers a computational speedup of order
$\tau_{p}/\tau_{0}$, which is large in extremely persistent active matter.
Second, the intermittent mechanism of ADD relaxation – which involves long
quiescent periods punctuated by large sudden motions – means that these
systems share similarities with other physical systems where athermal
quasistatic motion is relevant. These include amorphous solids under athermal
quasistatic shear (AQS) maloney2006amorphous , and active matter systems
undergoing an athermal quasistatic random displacement (AQRD) protocol
morse2021direct .

In all these contexts, a system evolves in response to potential gradients and
the potential evolves quasistatically: the resulting motion is mostly
quasistatic, but sudden relaxation events are triggered when the local minimum
of the potential deforms into a saddle point. Such events are known as
avalanches, because large-scale motion can be triggered after an infinitesimal
local change sethna1993hysteresis ; liu2016driving ; ozawa2018random .

In addition to similarities with AQS and AQRD, the ADD method is also related
to the athermal quasistatic random force (AQRF) protocol of Ref.
morse2021direct . However, ADD is distinct from all these methods.
Specifically, both AQS and AQRD displace particles along a fixed driving
direction and particles move orthogonally to the drive to minimise their
interaction energy. This leads to sustained stick-slip motion, including
avalanches. By contrast, AQRF applies forces with fixed direction; these are
changed smoothly, which generates avalanches. However, increasing the force
amplitude eventually drives the system through a yielding threshold
liao2018criticality ; mandal2020extreme ; villarroel2021critical , after which
mechanical equilibrium is no longer reached and the system “flows”. The
essential features of ADD are that it controls the forces on particles
(contrary to AQS and AQRD); and that the directions of these forces change
randomly with time while their typical strength remains constant, contrary to
AQRF. The result is a dynamical non-equilibrium steady state that exhibits
stick-slip motion. This corresponds to the large-$\tau_{p}$ limit of the
steady states previously observed in active matter, which also show highly
intermittent motion mandal2020extreme ; keta2022disordered .

In this work, we develop an efficient computational approach to exploit the
ADD method in order to characterise the dynamical relaxation of dense active
matter over timescales $t\gg\tau_{p}$ in the extremely persistent limit of
$\tau_{p}\gg\tau_{0}$. We work in two dimensions, which is the relevant case
for many experimental active systems angelini2011glasslike ; kumar2014flocking
; dauchot2016crystal , and also helps visualisation of the complex physics of
active systems. We address two broad sets of questions.

First, we show that ADD relaxation is built on two fundamental processes as in
other athermal quasistatic methods such as AQS. These are smooth elastic
deformations where the particles deform weakly near an energy minimum, in
response to the slowly-varying self-propulsion forces; and sudden plastic
rearrangements or avalanches in which a local minimum of the energy landscape
becomes unstable, forcing the system towards a new one. We characterise these
two processes in detail, including quantitative comparisons with AQS.

Our second broad question is how particles move and relax in the steady states
of ADD. As one might expect for a dense disordered system with slow dynamics,
one finds cooperative heterogeneous relaxation. We analyse this motion using a
range of techniques borrowed from passive glassy systems berthier2011dynamical
, including distributions of particle displacements, correlation functions
tailored for structural relaxation, and four-point susceptibilities.

Our results generate significant insight into active matter at large
$\tau_{p}$. For the elastic and plastic events, several of our results show
spatial correlations that extend across the entire system. For example, the
mean size of a plastic event grows with system size $N$ mandal2021how , and we
also observe a significant $N$-dependence of the yielding threshold
mandal2020extreme . This should be contrasted with thermal glassy systems
where one expects the behaviour of a large bulk system to be well-approximated
by dividing it into several non-interacting subsystems. This behaviour appears
because we choose to take the limit $\tau_{p}\to\infty$ at fixed system size,
which implies that the characteristic length scale for velocity correlations
grows with the system size henkes2020dense . We discuss the consequences of
this choice when comparing ADD dynamics with direct simulation of active
matter.

For structural relaxation mechanisms in the steady state, we find behaviour
with many similarities to passive glasses, including dynamical slowing down,
stretched exponential relaxation, and significant dynamical heterogeneity,
including dynamical facilitation effects chandler2011 . On the other hand, the
standard two-step relaxation mechanism of thermal glasses berthierbiroli11
and weakly persistent active glasses berthier2019glassy  is not observed.
There is no notion of thermal motion within a cage and particles instead move
slowly and smoothly in elastic trajectory segments before suddenly relaxing by
plastic ones. The second major distinction with thermal systems is the
existence of system-spanning dynamic correlations that control the behaviour
of time correlations functions in a way that is again qualitatively
reminiscent of AQS in sheared glasses.

The remainder of the paper is organised as follows. Sec. II defines the model
and presents our implementation of the ADD method. Sec. III characterises the
elastic and plastic steps. Sec. IV analyses the structural relaxation and Sec.
V collects our conclusions.

##  II Model and methods

###  II.1 Model

We consider a system of $N$ polydisperse overdamped active Ornstein-Uhlenbeck
particles (AOUPs) in a two-dimensional periodic box of linear size $L$. This
is the same system as in keta2022disordered  where large finite $\tau_{p}$’s
were considered; here we use ADD to access the limit $\tau_{p}\to\infty$.
Particle $i$ has position $\bm{r}_{i}$ and diameter $\sigma_{i}$; it feels a
self-propulsion force $\bm{p}_{i}$ with persistence time $\tau_{p}$; and it
interacts with other particles through a regularised soft-sphere Weeks-
Chandler-Andersen potential of strength $\varepsilon$ (see SM). The resulting
equations of motion are

$\displaystyle\xi\dot{\bm{r}}_{i}(t)$
$\displaystyle=-\nabla_{i}U(t)+\bm{p}_{i}(t),$ (1)
$\displaystyle\tau_{p}\dot{\bm{p}}_{i}(t)$
$\displaystyle=-\bm{p}_{i}(t)+\sqrt{2\xi^{2}D_{0}}\,\bm{\eta}_{i}(t),$ (2)

where $\xi$ is a viscous damping coefficient, $D_{0}$ is the free-particle
diffusion constant, and each component of the vector $\bm{\eta}_{i}$ an
independent zero-mean unit-variance Gaussian white noise. The propulsions
$\bm{p}_{i}$ follow a non-interacting Ornstein-Uhlenbeck process, and thus
their distribution is Gaussian at each time. Particle diameters are drawn from
a uniform distribution of mean $\sigma=\overline{\sigma}_{i}$ and
polydispersity 20% fily2014freezing ; keta2022disordered . The motivation for
studying polydisperse systems is to suppress the tendency of these systems to
crystallise digregorio2018full , and to allow access to the dense fluid regime
where motion is complex and cooperative keta2022disordered .

At fixed propulsion forces, the system converges to a force-balanced
configuration in a time $\approx\tau_{0}=\xi\sigma^{2}/\varepsilon$ if these
forces are not too strong mandal2021how . We focus here on extremely
persistent systems where $\tau_{p}\gg\tau_{0}$. In this limit, it is
convenient to introduce a rescaled time variable as

$t^{\prime}=t/\tau_{p}\;.$ (3)

We also rescale the propulsion by defining $\tilde{\bm{p}}=\bm{p}/f$ with
$f=\sqrt{\xi^{2}D_{0}/\tau_{p}}$. It is convenient to work in the centre-of-
mass frame so we write
$\overline{\tilde{\bm{p}}}=(1/N)\sum_{i}\tilde{\bm{p}}_{i}$ and
$\overline{\bm{r}}=(1/N)\sum_{i}\bm{r}_{i}$, which are the propulsive force
acting on the particles’ centre of mass and the position of the latter. In the
centre-of-mass frame and using the above reduced units, Eqs. (1,2) become

$\displaystyle\frac{\xi}{\tau_{p}}\frac{\mathrm{d}\bm{r}^{\prime}_{i}}{\mathrm{d}t^{\prime}}(t^{\prime})$
$\displaystyle=-\nabla_{i}U(t^{\prime})+f[\tilde{\bm{p}}_{i}(t^{\prime})-\overline{\tilde{\bm{p}}}(t^{\prime})],$
(4)
$\displaystyle\frac{\mathrm{d}\tilde{\bm{p}}_{i}}{\mathrm{d}t^{\prime}}(t^{\prime})$
$\displaystyle=-\tilde{\bm{p}}_{i}(t^{\prime})+\sqrt{2}\,\bm{\eta}^{\prime}_{i}(t^{\prime}),$
(5)

where $\bm{r}^{\prime}_{i}=\bm{r}_{i}-\overline{\bm{r}}$ is the position
relative to the centre of mass, and $\bm{\eta}^{\prime}_{i}$ is a zero-mean
Gaussian white noise in the rescaled time variables, that is,
$\bm{\eta}^{\prime}_{i}=\sqrt{\tau_{p}}\bm{\eta}_{i}$, which ensures that the
components of $\bm{\eta}^{\prime}_{i}$ satisfy
$\langle\eta^{\prime}_{i}(t_{1}^{\prime})\eta^{\prime}_{i}(t_{2}^{\prime})\rangle=\delta(t_{1}^{\prime}-t_{2}^{\prime})$.

To arrive at the ADD limit, we take limit $\tau_{p}\to\infty$ at fixed $f$.
For large $\tau_{p}$, Eq. (4) describes very fast relaxation to configurations
with perfect force balance, which satisfy

$\nabla_{i}U_{\rm eff}=0,$ (6)

where

$U_{\rm
eff}=U-f\sum_{j}[\tilde{\bm{p}}_{j}(t^{\prime})-\overline{\tilde{\bm{p}}}(t^{\prime})]\cdot\bm{r}_{j}$
(7)

is an effective potential, which corresponds to the original potential energy
of the interacting particle system, tilted by the active forces. Hence, as
$\tau_{p}\to\infty$, the system is almost always in a local minimum of $U_{\rm
eff}$, as captured mathematically by Eq. (6).

It is natural to fix $\sigma$ as the unit of length and $\varepsilon$ as the
unit of energy. The rescaled time $t^{\prime}$ is already dimensionless. In
the ADD limit $\tau_{p}\to\infty$ studied here, the only remaining control
parameters are therefore the number density $\rho=N/L^{2}$, the self-
propulsion force $f$, and the total number of particles $N$. We work
throughout at the representative density $\rho=1.2$. For fixed $f$, we expect
the behaviour to be robust with respect to $\rho$, as long as $\rho$ is not
too large (complete jamming) or too small (no force balanced states, i.e.
flow).

###  II.2 The quasi-static ADD method

The dynamics of our model can be simulated directly in the ADD limit
mandal2021how ; we summarise here the method for achieving this. The steady
state distribution of the propulsions $\tilde{\bm{p}}_{i}$ factorises across
particles, with

$P(\tilde{\bm{p}}_{i})=\frac{1}{2\pi}\exp\left(-\frac{1}{2}|\tilde{\bm{p}}_{i}|^{2}\right),$
(8)

and we use this to initialise the $\tilde{\bm{p}}_{i}$. In each step of the
ADD simulation, the propulsion dynamics in Eq. (5) is first integrated with
time step $\delta t^{\prime}$:

$\tilde{\bm{p}}_{i}(t^{\prime}+\delta t^{\prime})=(1-\delta
t^{\prime})\tilde{\bm{p}}_{i}(t^{\prime})+\sqrt{2\,\delta
t^{\prime}}\,\tilde{\bm{\eta}}^{\prime}_{i},$ (9)

where
$\tilde{\bm{\eta}}^{\prime}_{i}=(\tilde{\eta}^{\prime}_{i,x},\tilde{\eta}^{\prime}_{i,y})$
are two random numbers drawn from a Gaussian distribution with zero mean and
unit variance.

Next one integrates Eq. (4), which requires that $U_{\rm eff}$ is minimised by
steepest descent, holding the $\tilde{\bm{p}}_{i}$ fixed. To increase
computational efficiency, we replace this steepest descent by a faster
conjugate gradient minimisation, using the GPL-licensed ALGLIB C++ library
alglib . This is much more efficient than steepest descent, but it may
generally lead to different local minima nishikawa2022relaxation . We find
that the differences between conjugate gradient and steepest descent
algorithms are significant only in steps for which the system moves far from
its initial position. For these individual steps, we then automatically revert
to steepest descent, as originally proposed in Ref. mandal2021how . That is,
we choose a threshold in the mean-squared displacement
$(1/N)\sum_{i}|\Delta\bm{r}_{i}|^{2}>0.1$ to identify minimisation steps with
large total displacements. If this occurs during conjugate gradient
minimisation then we restart the minimisation step and use steepest descent
for that particular step. It is noteworthy that each of these minimisation
step requires on average $10^{2}$-$10^{3}$ force evaluations. It is thus
numerically challenging to explore both large systems and the large times
needed to reach steady state. We will therefore show data for $N\leq 2000$.

The ADD construction is valid for forces $f$ below an $N$-dependent yielding
threshold $f^{*}(N)$. The potential $U_{\rm eff}$ is not bounded below so
steepest descent may not converge to a local minimum – instead the particles
could continue to move along their self-propulsion directions. This happens
for $f>f^{*}$; the rheology and phase behaviour of the system above this
threshold is explored in Ref. villarroel2021critical . We locate the threshold
$f^{*}(N)$ liao2018criticality  by letting the system evolve from a random
initial configuration and checking for the average proportion of systems still
flowing after some time $t$. We obtain the following rough estimates:
$f^{*}(N=500)\approx 1.7$, $f^{*}(N=1024)\approx 1.4$, $f^{*}(N=2000)\approx
1.2$, $f^{*}(N=4096)\approx 1.1$. Given these systematic finite-size effects,
it would be desirable to simulate even larger systems. However, we emphasise
that a single time step for ADD is much more expensive than a single time step
in a standard molecular dynamics simulation, because it requires many
evaluations of the interparticle forces to converge the energy minimisation.
This stems from the intrinsic difficulty of simulating systems with well-
separated time scales ($\tau_{p}/\tau_{0}\to\infty)$. This limit is
inaccessible using standard methods: it can be simulated using ADD, but there
is still a significant cost.

In the ADD dynamics for $f<f^{*}$, each step starts with the system in a local
minimum of $U_{\rm eff}$ with positions and propulsions
$(\bm{r}^{0},\bm{p}^{0})$, and evolves to a new minimum with positions and
propulsions $(\bm{r},\bm{p})$. This can happen in two ways, as sketched in
Fig. 1(a). In the simplest case, a small change in propulsive forces changes
the local minimum of $U_{\rm eff}$ perturbatively, leading to small
displacements. This will be called an elastic step. However, the change in
propulsive forces can also destabilise the local minimum at
$(\bm{r}^{0},\bm{p}^{0})$, leading to a non-perturbative change in the
configuration. This is called a plastic step.

Figure 1: (a) Sketch of the effective potential energy landscape
$U_{\mathrm{eff}}$ in Eq. (7) at times 0 (thick blue line) and $\delta
t^{\prime}$ (thin red line). The system initially rests in a minimum of
$U_{\mathrm{eff}}(0)$ (filled blue circle). After $\delta t^{\prime}$, the
system rests in a minimum of the new landscape $U_{\mathrm{eff}}(\delta
t^{\prime})$, and it is displayed in the initial landscape (open red circle).
We distinguish elastic events ($\delta\epsilon_{p}\geq 0$) for which the
systems remains close to its original position, and thus the potential energy
in $U_{\mathrm{eff}}(0)$ increases, and plastic events
($\delta\epsilon_{p}<0$) for which a rearrangement occurs. (b) Accumulated
variations of the effective potential energy, $\sum_{t^{\prime\prime}\leq
t^{\prime}}\delta\epsilon_{p}(t^{\prime\prime})$. We identify elastic branches
made of successive elastic events (e.g. event 1) as ascending lines, and
plastic events (e.g. events 2 and 3) as instantaneous large drops. We define
$\tau^{\prime}$ as the time between two consecutive plastic events. Parameter
values: $f=0.9$, $N=1024$, $\delta t^{\prime}=10^{-2}$, purple square symbols
separated by $\delta t^{\prime}$.

To distinguish these two cases, we compute the change in effective potential
in one step:

$\delta\epsilon_{p}=U(\bm{r})-U(\bm{r}^{0})-f\sum_{i}(\tilde{\bm{p}}_{i}^{0}-\overline{\tilde{\bm{p}}^{0}})\cdot(\bm{r}_{i}-\bm{r}_{i}^{0}).$
(10)

This sign convention is opposite to that of Ref. mandal2021how . Note that the
propulsions in this equation are those of the state before the ADD step, with
the consequence that perturbative changes in the positions lead to positive
$\delta\epsilon_{p}$ as the system moves away from the minimum of the
associated $U_{\rm eff}$. Hence, we identify elastic steps as those with
$\delta\epsilon_{p}\geq 0$ while those with $\delta\epsilon_{p}<0$ are
plastic, as illustrated in Fig. 1(a). We discuss these two types of step
separately in detail below.

As well as the system parameters, a numerical simulation of ADD also requires
a choice of time step $\delta t^{\prime}$. As usual, this should be small
enough to mimic the limit $\delta t^{\prime}\to 0$, but large enough to ensure
computational efficiency. In practice, we chose $\delta t^{\prime}$ depending
on the state point: we tested several choices and monitored the mean squared
displacement during plastic events. We chose a value small enough that this
quantity depends at most weakly on $\delta t^{\prime}$ mandal2021how .

##  III Analysis of individual events

A typical trajectory from ADD is shown in Fig. 1(b). It consists of sequences
of elastic steps [$\delta\epsilon_{p}$ is positive and $O(\delta
t^{\prime})$], interspersed with instantaneous plastic events
[$\delta\epsilon_{p}$ is negative and $O(1)$]. Such behaviour is familiar from
AQS simulations of sheared glasses maloney2006amorphous  as well as from AQRD
simulations morse2021direct  and the non-equilibrium dynamics of the random-
field Ising model (RFIM, where the plastic events would be identified as
avalanches) ozawa2018random .

This section analyses the properties of the plastic and elastic steps,
including a comparison with AQS. We take $f=0.9$ throughout this section. This
is a practical choice: on the one hand it is far enough from the threshold
$f^{*}$ to keep the system from moving too much between minimisations; on the
other hand smaller values of $f$ lead to slower dynamics and the numerics
become more challenging, as discussed in Sec. IV.

###  III.1 Elastic steps

Figure 2: (a) Snapshot of displacements for a single elastic step (1 in Fig.
1). Colours indicate the norm and arrows their direction. Parameter values:
$f=0.9$, $N=1024$, $\delta t^{\prime}=10^{-2}$. (b) Corresponding displacement
correlation function
$C(r)=\left<\delta\bm{r}_{i}\cdot\delta\bm{r}_{j}\delta(r-r_{ij})\right>/\left<|\delta\bm{r}_{i}|^{2}\right>$
for various $N$ values. The collapse with $r/L$ for three different system
sizes shows that correlations scale with system size, becoming negative for
large enough $r/L$. Solid (resp. dotted) line: evaluation of (16) (resp. (17))
with sums across the range $0<m^{2}+n^{2}<40^{2}$ and a suitably chosen
overall prefactor.

A representative snapshot of the displacement field obtained during an
elementary elastic step of the ADD dynamics is shown in Fig. 2(a). We observe
highly heterogeneous displacements, with wide variations in amplitude and
clear large-scale correlations resembling both non-affine displacement in
sheared athermal glasses and collective swirling motion in active matter
henkes2020dense .

Elastic steps can be analysed under the assumption that the updated propulsive
forces move the minimum of $U_{\rm eff}$ perturbatively, in a way similar to
AQS maloney2006amorphous ; maloney2006correlations . The states before and
after the elastic step are both force balanced, so for all $i$

$\delta\left[-\nabla_{i}U+f(\tilde{\bm{p}}_{i}-\overline{\tilde{\bm{p}}})\right]=0,$
(11)

where $\delta$ indicates the change in a time increment $\delta t^{\prime}$.
For small $\delta t^{\prime}$ this implies that

$-\sum_{j}(\mathbb{H})_{ij}\delta\bm{r}_{j}+\bm{\Xi}_{i}=0,$ (12)

where $\mathbb{H}$ is the Hessian matrix of $U$, that is
$(\mathbb{H})_{ij}=\frac{\partial^{2}U}{\partial\bm{r}_{i}\partial\bm{r}_{j}}$,
and

$\displaystyle\bm{\Xi}_{i}$
$\displaystyle=f(\delta\tilde{\bm{p}}_{i}-\delta\overline{\tilde{\bm{p}}})$
(13) $\displaystyle=f\,(-\tilde{\bm{p}}_{i}\delta t^{\prime}+\sqrt{2\,\delta
t^{\prime}}\bm{\eta}_{i}^{\prime}-\delta\overline{\tilde{\bm{p}}})$

is the analogue of the affine force in the AQS setting of Ref.
maloney2006amorphous , where the term $\delta\overline{\tilde{\bm{p}}}$ only
ensures $\sum_{i}\bm{\Xi}_{i}=0$ so that we stay in the centre-of-mass frame.

To evaluate the solution of (12) approximately, we follow Ref.
maloney2006amorphous  in assuming that the eigenmodes of the Hessian
$\mathbb{H}$ can be decomposed as plane wave eigenmodes of the Navier operator
slaughter2002linearized  (see SM)

$\displaystyle\delta\bm{r}_{i}=$
$\displaystyle\sum_{j}(\mathbb{H}^{-1})_{ij}\bm{\Xi}_{j}$ (14)
$\displaystyle=$
$\displaystyle\sum_{m,n,\alpha}\frac{\Xi_{mn}^{\alpha}}{\lambda^{\alpha}(m^{2}+n^{2})}e^{\mathrm{i}\bm{k}_{mn}\cdot\bm{r}_{i}}\hat{\bm{k}}^{\alpha}_{mn},$

where $\Xi_{mn}^{\alpha}$ is the projection of the affine force along the
corresponding eigenmode, $\bm{k}_{mn}=(2\pi m/L,2\pi n/L)$ is the wave vector
and the third sum is over the two polarization directions of the elastic
displacements $\alpha=||,\perp$ (longitudinal and transverse, respectively).
The corresponding polarization vectors $\hat{\bm{k}}_{mn}^{\alpha}$ are the
unit vectors parallel and orthogonal to the wave vectors, respectively. We
have written the associated eigenvalues as $\lambda^{\alpha}(m^{2}+n^{2})$,
with the prefactor again depending on the polarization direction. It is then
possible to compute the spatial correlation function

$\displaystyle g^{\mathrm{ADD}}(\bm{r})\propto$
$\displaystyle\left<\delta\bm{r}_{i}\cdot\delta\bm{r}_{j}\,\delta(\bm{r}-(\bm{r}_{j}-\bm{r}_{i}))\right>$
(15) $\displaystyle=$
$\displaystyle\sum_{m,n,\alpha}\frac{\left<|\Xi_{mn}^{\alpha}|^{2}\right>}{(\lambda^{\alpha}(m^{2}+n^{2}))^{2}}e^{\mathrm{i}\bm{k}_{mn}\cdot\bm{r}},$

where we have checked numerically that the projections of the affine force
$\bm{\Xi}_{i}$ on the eigenmodes of the Hessian $\mathbb{H}$ behave as
uncorrelated random numbers, with a variance
$\left<|\Xi_{mn}^{\alpha}|^{2}\right>$ independent of the specific mode. We
finally take the orientational average of (15) to write

$g^{\mathrm{ADD}}(r)\propto\sum_{m,n}\frac{J_{0}\left(2\pi\sqrt{m^{2}+n^{2}}\,r/L\right)}{(m^{2}+n^{2})^{2}},$
(16)

where $J_{0}$ is the zeroth-order Bessel function of the first kind. Our
result is quantitatively different from the corresponding correlation function
for sheared amorphous solids maloney2006correlations

$g^{\mathrm{AQS}}(r)\propto\sum_{m,n}\frac{J_{0}\left(2\pi\sqrt{m^{2}+n^{2}}\,r/L\right)}{m^{2}+n^{2}}\;.$
(17)

The difference between (16,17) arises because affine forces $\bm{\Xi}_{i}$ in
AQS are derived from pair potentials: the force exerted by particle $i$ on
particle $j$ is equal and opposite to the one exerted by particle $j$ on
particle $i$. Hence, the affine forces on different particles are necessarily
correlated in AQS. Still, both functions correspond to scale-free
correlations, with the only relevant length scale being the system size $L$
itself. Qualitatively, this behaviour is visible already from the system-
spanning vortices in the displacement field shown in Fig. 2(a).

We compute the displacement correlation quantitatively in Fig. 2(b). These
correlations scale with the system size for the range of sizes we have
investigated, and are close to the analytical prediction (16). The small
difference between the prediction and the measurement may be attributed to the
plane wave hypothesis, which can in principle be tested mizuno2017continuum .
We have also checked that the above correlation functions are consistent with
the correlations computed from purely harmonic steps, i.e. with displacements
determined by solving (12) exactly. This establishes another parallel to Ref.
maloney2006correlations .

The scaling with system size of the displacement correlations that we find is
consistent with the arguments obtained for finite persistence time in Ref.
henkes2020dense . In that case, the dynamics along elastic trajectory segments
produces displacement (or equivalently velocity) correlations on a length
scale that diverges as $\sim\sqrt{\tau_{p}}$ for large persistence times. As
our analysis takes $\tau_{p}\to\infty$ from the start, this limit translates
into displacement correlations on the largest length scale available, i.e. the
system size.

###  III.2 Plastic steps

Figure 3: (a,b) Snapshots showing the movement of particles in the two plastic
events marked 2 and 3 in Fig. 1, with participation $S=19$ and $S=102$,
respectively. Displacements are magnified 5 times and superimposed onto a
colormap of the residual force $|\bm{f}_{\mathrm{res},i}|$ that highlights
rearranging particles. Parameter values: $f=0.9$, $N=1024$, $\delta
t^{\prime}=10^{-2}$. (c) Log-distribution of the participation $S$ in plastic
events, for three system sizes $N$ and time steps $\delta
t^{\prime}(N=500)=2\times 10^{-3}$, $\delta t^{\prime}(N=1024)=10^{-3}$,
$\delta t^{\prime}(N=2000)=5\times 10^{-4}$. The solid black line corresponds
to a scaling $P(S)\sim P(\log_{10}S)/S\sim S^{0.3-1}=S^{-0.7}$. Inset shows
the evolution of the mean $\langle S\rangle$ with $N$. (d) Log-distribution of
times $\tau^{\prime}$ between consecutive plastic events. Inset shows the
evolution of the mean $\langle\tau^{\prime}\rangle$ with $N$. Scaling
exponents are subject to significant uncertainties and the numbers provided
are indicative.

Particle displacements for two representative plastic steps are shown in Figs.
3(a,b). The qualitative picture has again many similarities with AQS: the
displacements in any single plastic event can be interpreted as sequences of
localised yielding events mandal2021how ; maloney2006amorphous . That is, a
plastic step happens when a local minimum of $U_{\rm eff}$ develops an
unstable direction, causing local motion; but the elastic perturbation due to
this event perturbs the system over large length scales and can create further
unstable directions in other parts of the system. This leads to a cascade or
avalanche of localised yielding events. We describe these plastic events here
and discuss similarities and differences with AQS.

We find that plastic avalanches display a broad range of sizes and can involve
a few localised particles, as in Fig. 3(a), or a greater number of particles
distributed across the system as in Fig. 3(b), or even the whole system. To
characterise the participation in each plastic event, we identify the number
$S$ of particles with significant changes in their local environments. In Ref.
mandal2021how  this identification was carried out by thresholding particle
displacements, but such a criterion neglects the fact that particles may
collectively move large distances without changing their local environment.
Here, we use instead the residual force $\bm{f}_{\mathrm{res},i}$
lemaitre2014structural ; lerbinger2021relevance  (see SM for definition) as an
indicator of rearrangements. This force is zero if and only if displacements
result from a harmonic response of the system to the change in propulsion
forces. If a localised avalanche takes place, particles that are far from the
avalanche tend to respond elastically, leading to very low residual forces,
whereas rearranging particles in the core of the avalanche have large residual
forces. This is illustrated in Fig. 3(a,b). We then define the avalanche size
$S$ as the number of particles for which $|\bm{f}_{\mathrm{res},i}(\delta
t^{\prime})|>20$, determined after careful analysis of the distributions of
$\bm{f}_{\mathrm{res},i}$ (see SM).

Fig. 3(c) shows the resulting broad distribution of $\log_{10}S$ for three
system sizes. The small events can be attributed to local yielding, in which
the remainder of the system reacts elastically with particles moving
collectively to accommodate the local rearrangement mandal2021how ;
lerbinger2021relevance . At low values of $S$, the distributions overlap for
different system sizes, with a behaviour compatible with $P(S)\sim S^{-\tau}$
with $\tau\approx 0.7$. (We use the notation $\tau$ here, which is standard in
the literature karimi2017inertia ; $\tau$ does not indicate a time scale.)
Given the small range of system sizes studied here, it is difficult to provide
a very precise estimate of $\tau$, but it is clearly distinct from the values
found in AQS simulations of sheared glasses where values in the range
$\tau\approx 1.2-1.5$ have been reported liu2016driving ; karimi2017inertia .

Turning to the behaviour at large $S$ we observe that larger avalanches with
$S\sim N$ are more frequent for larger systems, suggesting that these are also
important in ADD. Correspondingly, the inset in Fig. 3(c) shows that the
average event size $\langle S\rangle$ scales as $N^{\gamma}$ with
$\gamma\approx 0.7$, showing that the mean avalanche size is indeed controlled
by large avalanches that are limited by the system size only. In other words,
the avalanches observed during plastic events also lead to scale-free dynamic
relaxation events.

Fig. 3(d) shows the distribution of waiting times $\tau^{\prime}$ between
consecutive events. The average time decreases with system size as
$\langle\tau^{\prime}\rangle\sim N^{-1.5}$. If localised yielding events
happened independently in different parts of a large system, one would have a
more trivial dependence on the system size, $\langle\tau^{\prime}\rangle\sim
N^{-1}$: together with the $N$-dependence of $\langle S\rangle$, this is
another indication of long-ranged correlations, on the scale of the system
size.

Such scaling behaviour hints at critical phenomena. The force threshold
$f^{*}(N)$ for yielding decreases with $N$ in our simulations, which is
presumably also due to long-range correlations liao2018criticality ;
mandal2020extreme . Since we increase $N$ at fixed $f$ while staying always
below $f^{*}(N)$, some of the dependence on $N$ may arise because the larger
systems are closer to yielding. Indeed, larger systems support larger events,
which tend to relax the system more quickly; compare Fig. 7(b) below.

While the ADD plastic events share similarities with those of AQS, there are
also some important differences to emphasise. In particular, in ADD there is
no preferred direction and the system is isotropic (apart from a possible
influence of the periodic boundary conditions, which we expect to be very
weak). In AQS, on the other hand, rotational symmetry is broken because the
system is always sheared in the same direction. As a result, localised plastic
events eventually organise into a line of slip, which leads to a subextensive
scaling of event sizes $\left<S\right>\sim L\sim\sqrt{N}$ in the steady state
maloney2004subextensive ; maloney2006amorphous . These correlations also cause
a reduction in the frequency of plastic events: the typical time (accumulated
strain) between consecutive events scales as $\left<\tau^{\prime}\right>\sim
1/L\sim 1/\sqrt{N}$ (we recall it would be $\sim 1/N$ for independent events
maloney2006amorphous ).

In short, ADD in its elastic steps produces displacements that are correlated
on the scale of the system size $L$, in agreement with predictions for finite
$\tau_{p}$ henkes2020dense . However, the master curve for displacement
correlations against $r/L$ is different from the AQS case
maloney2006correlations  because the local “affine” forces $\bm{\Xi}_{i}$ lack
the correlations that are present for AQS maloneyprivate . For plastic events,
the isotropy of ADD also leads to larger event sizes $S$ and shorter inter-
event times. The next question to be addressed is how the individual particles
move in the active fluid, when observed over multiple time steps.

##  IV Microscopic dynamics

Since ADD is a computational tool to explore particle motion in dense active
fluids, it is natural to study dynamical relaxation in ADD trajectories. To
this end, we use observables developed for the analysis of relaxation in dense
glassy systems berthier2011dynamical . Such measurements have also been used
to describe particle motion in sheared and active glasses, which all display
heterogeneous and cooperative dynamics.

###  IV.1 Mean squared displacement

Figure 4: (a) Mean-squared displacement for different values of the self-
propulsion force $f$. (b) Distribution of displacements scaled by the MSD at
different $t^{\prime}$. Parameter values: $N=1024$, $\delta
t^{\prime}=10^{-2}$, $f=0.9$.

Fig. 4(a) shows the mean squared displacement (MSD)

$\mathrm{MSD}(t^{\prime})=\left<|\bm{r}_{i}(t^{\prime})-\bm{r}_{i}(0)|^{2}\right>$
(18)

for different values of the self-propulsion force $f$. In the steady state,
the MSD is nearly diffusive at all times. The self-diffusion constant at large
times roughly decreases by an order of magnitude between $f=0.9$ and $f=0.8$.
However, there is no feature in the average displacements that would allow
identification of a characteristic relaxation time scale or length scale. This
is in contrast to the classic two-step relaxation scenario found in many
glassy systems keta2022disordered ; binder2005glassy ; berthierbiroli11 , but
resembles the diffusive behaviour found in AQS simulations of sheared systems
heussinger2010superdiffusive .

Figure 5: (a) Elastic and plastic contributions to the MSD, together with the
MSD predicted from a Gaussian fit of the central part of $G_{s}$. (b) Wave-
vector-dependent relaxation time $\tau^{\prime}_{\mathrm{s}}(k)$ extracted
from elastic and plastic displacements (colours and symbols as in (a)).
Parameter values: $f=0.9$, $N=1024$, $\delta t^{\prime}=10^{-2}$.

Although the MSD displays seemingly trivial behaviour, the displacement
distributions have significant structure. Fig. 4(b) shows the corresponding
distribution $G_{s}(r)$ of the $x$ and $y$-components of the particle
displacements, scaled by the root mean-squared displacement, at $f=0.9$. These
distributions differ strongly from Gaussian behaviour, which is only recovered
in the large time limit, $t^{\prime}\to\infty$. At small times $t^{\prime}\ll
1$ (with of course $t^{\prime}\geq\delta t^{\prime}$), the displacement
distribution has a narrow central peak with heavy tails. The width of these
tails decreases with increasing time, and the distribution approaches a
Gaussian form. For supercooled liquids, we would expect the small time
distribution to be nearly Gaussian due to short-time thermal dynamics, with
fat tails developing only as the system starts to relax. The tails appear when
a significant number of particle rearrangements has taken place
kob1997dynamical ; chaudhuri2007universal . The difference between ADD and
thermal dynamics at short times is easily explained by the athermal
quasistatic nature of ADD dynamics. Moreover, the participation in plastic
events has a broad distribution [Fig. 3(c)]. As a result, the fat tails
arising from structural relaxation are visible already after a single time
step $\delta t^{\prime}$ and arise from avalanches. In AQS simulations,
similar heavy-tailed distributions also appear at early times due to plastic
avalanches heussinger2010superdiffusive . In both ADD and AQS we expect that
nearly-exponential tails appear at intermediate times, as a generic result of
the stochastic nature of avalanches chaudhuri2007universal .

To gain more insight into the dynamics, we decompose the displacements into
separate contributions from elastic and plastic events. We define the elastic
(resp. plastic) displacement of a particle between $0$ and $t^{\prime}$ as the
sum of its displacements over all elastic (resp. plastic) steps between these
two times. We plot in Fig. 5(a) the MSDs from these contributions at $f=0.9$.
Both of them show a crossover between two diffusive scaling regimes, i.e. both
have ${\rm MSD}(t^{\prime})\sim t^{\prime}$ at short and at long times but
with different prefactors. Despite the complex time dependences of the
separate contributions, their sum in the total MSD appears nearly linear
(recall Fig. 4).

To connect the elastic and plastic displacements to the distribution
${G}_{s}(r)$ in Fig. 4(b), we fit the central peak of ${G}_{s}(r)$ to a
Gaussian distribution with standard deviation $\mathrm{std}(t^{\prime})$ such
that the associated mean-squared displacement is
$\mathrm{MSD}(t^{\prime})=2\,\mathrm{std}^{2}(t^{\prime})$. Fig. 5(a) compares
this effective MSD to the elastic and plastic contributions. At small times,
the central peak of ${G}_{s}$ is compatible with the variance of elastic
displacements, while at large times it is compatible with plastic
displacements. Our interpretation is that for small times, displacements in
elastic branches populate the narrow central part of the displacement
distribution, while the displacements of rearranging particles populate the
tails. At large times, plastic displacements dominate elastic displacements,
and the whole displacement distribution is close to Gaussian, therefore its
variance is dictated by plastic displacements.

###  IV.2 Single-particle correlation functions

The MSD and the distribution $G_{s}$ yield useful information about the
dynamics. However, MSD measurements can sometimes be dominated by a subset of
fast moving particles. To investigate this, we computed the self-intermediate
scattering function binder2005glassy

$F_{s}(k,t^{\prime})=\big{<}\cos\big{[}\bm{k}\cdot(\bm{r}_{i}(t^{\prime})-\bm{r}_{i}(0))\big{]}\big{>}$
(19)

where $\bm{k}$ is a suitable wavevector and $k=|\bm{k}|$. We analyse the
$k$-dependent relaxation time scale $\tau^{\prime}_{s}(k)$ berthier2005length
; berthier2004time , which is defined as the time at which
$F_{s}(k,\tau^{\prime}_{s}(k))=1/{\rm e}$.

We plot $\tau^{\prime}_{s}(k)$ for both elastic and plastic displacements in
Fig. 5(b). Fickian diffusion would correspond to the scaling
$\tau^{\prime}_{s}(k)\sim k^{-2}$. Elastic displacements show two distinct
diffusive scalings, at small and at large length scales (resp. small and large
time scales). At small times, elastic displacements are computed over a single
elastic trajectory segment, thus the corresponding diffusive behaviour derives
from the balance between the Ornstein-Uhlenbeck driving on the one hand and
the restoring forces on the other hand henkes2012extracting ;
bottinelli2017how . At large times, elastic displacements are computed over
many elastic trajectory segments separated by multiple plastic events. We
expect the displacements over these different elastic trajectory segments to
be independent, therefore the sum of all these displacements produces a
diffusive behaviour distinct from that of single branch displacements.

The small and large length scale behaviour are also different for plastic
displacements. These displacements are indeed diffusive at large length scales
but they show a relaxation time scale $\tau^{\prime}_{s}(k)$ that plateaus at
small length scales. This plateau corresponds to the typical time for a
plastic event to occur, and can be computed from the statistics of the inter-
event times $\tau^{\prime}$ as the residual time
$\tau^{\prime}_{\mathrm{res}}=\langle{\tau^{\prime}}^{2}\rangle/(2\langle\tau^{\prime}\rangle)$
pal2022inspection . Therefore, at times
$t^{\prime}\lesssim\tau^{\prime}_{\mathrm{res}}$ the plastic MSD is likely
dominated by a small subset of particles.

It is noteworthy that this result is robust to changes in $f$ (data not
shown). The residual time $\tau^{\prime}_{\mathrm{res}}$ changes by a factor
of $\sim 2.5$ between $f=0.9$ and $f=0.8$, distinct from the factor of $10$
observed for the self-diffusion constant (not shown). Moreover, since the
diffusion constant drops more rapidly than the residual time increases, the
typical length scale above which the plastic movement appears Fickian
berthier2005length  decreases with decreasing $f$ – which is opposite to what
we would expect for a supercooled liquid approaching the glass transition.

Figure 6: Bond breaking correlation function $C_{\mathrm{b}}(t^{\prime})$ for
different values of the self-propulsion force $f$. Parameter values: $N=1024$,
$\delta t^{\prime}=10^{-2}$.

The displacement fields in the plastic events of Figs. 3(a,b) show that
particles that are not involved in rearrangements can move away from their
initial position without relaxing their local structure. But neither the MSD
nor $F_{s}(k,t^{\prime})$ can detect whether single particle translations
actually correspond to changes in the local structure or not. To focus on this
aspect of structural relaxation, we use the bond breaking correlation function
$C_{\mathrm{b}}(t)$ shiba2012relationship . Denoting
$\hat{r}_{ij}(t^{\prime})=|\bm{r}_{j}(t^{\prime})-\bm{r}_{i}(t^{\prime})|/\sigma_{ij}$
as the rescaled distance between particles $i$ and $j$ at time $t^{\prime}$,
with $\sigma_{ij}=(\sigma_{i}+\sigma_{j})/2$, we define

$C_{\mathrm{b}}(t^{\prime})=\frac{\sum_{i,j}\Theta(A_{1}-\hat{r}_{ij}(0))\,\Theta(A_{2}-\hat{r}_{ij}(t^{\prime}))}{\sum_{i,j}\Theta(A_{1}-\hat{r}_{ij}(0))}\;,$
(20)

where the parameter $A_{1}=1.25$ is a cutoff defining initial neighbours,
$A_{2}=1.5$ quantifies the distance they are required to separate before the
correlation function decays, and $\Theta$ designates the Heaviside function.
This function obeys $C_{b}(t^{\prime}=0)=1$, by definition, and it quantifies
at time $t^{\prime}$ the average fraction of neighbours lost since
$t^{\prime}=0$. This way, it efficiently disentangles rearrangements from
displacements that do not relax the local structure.

Fig. 6 shows $C_{\mathrm{b}}(t^{\prime})$ for several values of $f$. The
relaxation time scale $\tau_{b}^{\prime}$ of $C_{\mathrm{b}}$ roughly
increases by an order of magnitude between $f=0.9$ and $f=0.8$, mirroring the
decrease of the self-diffusion constant. Moreover, $\tau_{b}^{\prime}\gg 1$,
so structural relaxation happens long after self-propulsion forces have fully
decorrelated from their initial values which occurs for $t^{\prime}\sim 1$.
The correlation function is stretched, which suggests that structural
relaxation is temporally heterogeneous binder2005glassy . In addition, it is
remarkable that the MSD at the time $\tau_{b}^{\prime}$ where local structure
becomes fully decorrelated is greater than unity. This is again very different
from thermal glasses where the escape from the cage also coincides with
structural relaxation. Here instead particles travel comparatively larger
distance without necessarily relaxing the structure, which can be seen as a
consequence of the swirling motion observed in snapshots such as Fig. 2(a)
(for elastic events) or Fig. 3(a,b) (for plastic ones).

Figure 7: (a) Mean-squared displacement and (b) bond-breaking correlation
function for different $N$. Parameter values: $f=0.9$, $\delta
t^{\prime}(N=500)=2\times 10^{-3}$, $\delta t^{\prime}(N=1024)=1\times
10^{-3}$, $\delta t^{\prime}(N=2000)=5\times 10^{-4}$.

We finally show the dependence of the dynamics on system size $N$. Fig. 7
shows the MSD and the bond breaking correlation function for a fixed $f$ but
different values of $N$. For all values of $N$ studied, the MSD is diffusive
at large times and the relaxation of $C_{\mathrm{b}}(t^{\prime})$ is
stretched. Strikingly, however, both functions strongly depend on the system
size with no sign of a saturation at some large $N$ value. As $N$ increases,
the particles move faster, resulting in an MSD that is larger and a time
correlation function that decreases faster. Since we have simulated three
system sizes, it is not easy to determine a precise scaling of the self-
diffusion constant and the relaxation time $\tau^{\prime}_{b}$ with $N$. There
are two sources for the system size dependence of the dynamics as discussed
above: the frequency and size of the plastic events both increase with $N$,
recall Fig. 3. Thus, in the ADD regime, the dynamics is always sensitive to
the system size, as a result of the large persistence time limit. This is
again in good analogy with the AQS dynamics where the self-diffusion constant
also changes with system size lemaitre2007plastic .

Figure 8: (a) Bond breaking dynamical susceptibility $\chi_{b}(t^{\prime})$
for different self-propulsion forces $f$ with $N=1024$ and $\delta
t^{\prime}=10^{-2}$. (b) Bond breaking dynamical susceptibility
$\chi_{b}(t^{\prime})$ for different $N$ with $\delta
t^{\prime}(N=500)=2\times 10^{-3}$, $\delta t^{\prime}(N=1024)=1\times
10^{-3}$, $\delta t^{\prime}(N=2000)=5\times 10^{-4}$. (c-f) Snapshots of the
system highlighting the local bond breaking correlation $C_{\mathrm{b},i}$
between time $t^{\prime}=0$ and a lag time of (c) $t^{\prime}=3.86$ (124
plastic events, $C_{\mathrm{b}}=0.89$) (d) $t^{\prime}=5.02$ (196 plastic
events, $C_{\mathrm{b}}=0.65$) (e) $t^{\prime}=6.53$ (252 plastic events,
$C_{\mathrm{b}}=0.52$) (f) $t^{\prime}=8.49$ (355 plastic events,
$C_{\mathrm{b}}=0.41$). Parameter values: $N=2000$, $f=0.9$, $\delta
t^{\prime}=5\times 10^{-4}$.

###  IV.3 Dynamical heterogeneity

In glassy fluids, one generally expects complex heterogeneous dynamics, where
spatial fluctuations around the average dynamical behaviour are important for
understanding the relaxation dynamics berthier2011dynamical . The system
considered here also has this feature. It is illustrated in the snapshots of
Figs. 8(c-f) which show maps of a single-particle analogue of the bond-
breaking correlation function (20). This is defined as

$C_{\mathrm{b},i}(t^{\prime})=\frac{\sum_{j}\Theta(A_{1}-\hat{r}_{ij}(0))\,\Theta(A_{2}-\hat{r}_{ij}(t^{\prime}))}{\sum_{j}\Theta(A_{1}-\hat{r}_{ij}(0))},$
(21)

and represents the fraction of bonds of particle $i$ that have been broken up
to time $t^{\prime}$. Just like its global analogue, this function decays from
unity to zero as the environment of particle $i$ decorrelates from its initial
state.

The definition of a local relaxation function in Eq. (21) allows us to
visualise in real space how the initial structure of the system relaxes as
dynamics proceeds. The main observation in the snapshots of Fig. 8 is that the
spatial distribution of relaxed particles, at any given time, reveals strong
spatial correlations in the local dynamics.

While similar observations of spatially correlated dynamics are quite generic
in dense amorphous materials, the time series shown in Fig. 8(c-f) reveals
additional features beyond the mere existence of correlations. We observe that
at short times [Fig. 8(c)], only a few particles have relaxed, and the spatial
structure of the $C_{b,i}$ reveals the very few underlying plastic events that
have taken place in this particular trajectory. At larger times [Figs.
8(d,e)], one sees that additional structural relaxation events tend to happen
in the close vicinity of previous ones. These observations have also been made
in slowly relaxing supercooled liquids at equilibrium and this effect is known
as dynamical facilitation guiselin2022microscopic ; chandler2011 . At the
microscopic level, these effects must correspond to spatio-temporal
correlations between successive plastic events. Close to the relaxation time
[Fig. 8(e-f)], we can clearly identify fast regions where particles have been
involved in numerous rearrangements and have relaxed their local initial
structure. We also see slow regions where particles’ local environments remain
the same. Finally, at long times all particles have $C_{b,i}\approx 0$ and one
recovers a homogeneous picture. It is noteworthy that despite the absence of
shear to organise plastic activity into anistropic structures (namely, shear
bands maloney2006amorphous ), correlations spontaneously emerge between the
plastic centres.

The extent of these dynamical spatial correlations can be quantified via the
dynamical susceptibility shiba2012relationship ; berthier2011dynamical ;
lacevic2003spatially

$\chi_{\mathrm{b}}(t^{\prime})=N[\left<C_{b}(t^{\prime})^{2}\right>-\left<C_{b}(t^{\prime})\right>^{2}]\;.$
(22)

This function is plotted for different values of the self-propulsion force $f$
in Fig. 8(a), and for different numbers of particles $N$ in Fig. 8(b). As
usual, $\chi_{\mathrm{b}}$ shows a peak at the time close (but not exactly
equal) to ${\tau_{\mathrm{b}}^{\prime}}$ where $C_{b}$ starts to decrease
shiba2012relationship , indicating that the dynamics is most heterogeneous
around these times. The slowdown of the dynamics with decreasing $f$ or $N$
(Figs. 6, 7) is reflected in the corresponding increase of this peak time.
Moreover, the height of the peak tells us about the typical number of
particles involved in correlated clusters in Fig. 8 berthier2011dynamical . At
fixed $N$, the dynamical slowdown is accompanied by an increased cooperativity
of the relaxation, as is observed for liquids approaching the glass transition
lacevic2003spatially . This situation is different when the self-propulsion
force $f$ is kept fixed: the dynamics speeds up with increasing $N$ but it
also becomes more cooperative with a larger dynamical susceptibility. Recall
that ADD is defined by taking the limit $\tau_{p}\to\infty$ at fixed $N$: as a
consequence, the length scale that characterises velocity fluctuations is
slaved to the system size, and diverges for large $N$. The global correlation
function $C_{\mathrm{b}}$ and its fluctuations $\chi_{\mathrm{b}}$ both change
systematically with $N$, revealing that the long-time relaxation dynamics is
also sensitive to the system size, presumably because of a cooperativity
length scale that diverges with $N$.

These observations provide further evidence that spatially heterogeneous
dynamics is very generic in dense and disordered fluids. A major difference
with equilibrium supercooled liquids is the system size dependence observed
for the dynamical heterogeneity, indicating a diverging correlation length.
This is attributed to the quasi-static nature of the dynamics, as also found
in AQS simulations heussinger2010superdiffusive . In addition, the slow growth
of dynamic correlations with time (Fig. 8) reveals the role of dynamic
facilitation. Whereas facilitation has been described before in equilibrium
dynamics guiselin2022microscopic ; scalliet2022thirty ; PhysRevX.1.021013 ,
much less is understood about its consequences for sheared and active systems.
Our findings suggest that facilitation could also be a very generic feature
characterising the relaxation dynamics of dense and disordered fluids, and
this clearly deserves further study in the context of driven amorphous
materials.

##  V Conclusion

Although models of spherical self-propelled particles are among the simplest
models of active matter they exhibit rich physics with emerging structures,
phases and dynamical behaviours. A remarkable feature is the emergence of non-
trivial velocity correlations in fluid szamel2021long , glassy
szamel2015glassy ; henkes2020dense , and crystalline lorenzo20  states which
may, for non-arrested states, give rise also to interesting correlations of
the displacements and structural relaxation events.

In this context, the efficient implementation of activity-driven dynamics
(ADD) mandal2021how ; mandal2020multiple  enables us to study the relaxation
of dense systems of self-propelled particles in which the persistence time
$\tau_{p}$ is large compared to the microscopic time $\tau_{0}$ that the
system needs to reach an arrested state for a given set of self-propulsion
forces. On time scales $t=t^{\prime}\tau_{p}$ of the order of the persistence
time, the dynamics then becomes intermittent (Fig. 1). In the absence of
rearrangements, the system reacts elastically to changes in the self-
propulsion forces. The resulting movements are correlated on the length scale
of the system (Fig. 2). Consecutive elastic events may be interrupted by
plastic events that trigger instantaneous rearrangements. The participation in
these events has a broad distribution; outside of the plastic core forming
these avalanches, the remainder of the system moves collectively to
accommodate the rearranging regions (Fig. 3).

Relaxation of the whole structure happens through the accumulation of many of
these plastic events. This relaxation dynamics is nearly diffusive at all
times (Fig. 4) and spatially heterogeneous (Fig. 8), implying that plastic
events are not independent and tend to concentrate where they have already
happened, in a fashion reminiscent of dynamic facilitation.

Dense assemblies of self-propelled particles in two dimensions may serve as a
proxy for confluent cell tissues angelini2011glasslike  in the study of their
collective motion henkes2020dense . We also expect our results to be
transferable from two dimensions to three dimensions, as are the salient
features of the physics of glasses illing2017mermin  and cooperative motion in
dense active matter wysocki2014cooperative .

The limit of large persistence $\tau_{p}\to\infty$ is taken at fixed number of
particles $N$ and there is thus a dynamical length scale that scales with the
system size. As a consequence, the average dynamics of the system and its
fluctuations all depend on the system size (Figs. 7, 8).

Our computational abilities do not yet enable us to explore the $f\to 0$ limit
where the relaxation time would become large even in rescaled units
$t/\tau_{p}$. At fixed $N$, we expect that the typical time between plastic
events $\tau_{\mathrm{res}}$ would become $\gg\tau_{p}$ in the limit of small
$f$. We would then presumably recover a two-step relaxation scenario, with
$\beta$-relaxation corresponding to the diffusive elastic exploration of a
potential energy minimum, where the localisation length scale may itself
depend on the system size. Further explorations are needed to validate this
scenario and the nature of the very slow dynamics emerging in this limit. At
fixed $f$ we would also expect quantitative changes as $\rho$ is increased:
avalanches should become rarer but the evolution of their distribution remains
uncertain. Also interesting to study in future will be the connection between
plasticity in ADD and approaches to yielding in passive materials that argue
in favour of a mechanism based on fluctuating energy barriers NicMarBar14
rather than effective thermal activation over fixed barriers SolLeqHebCat97 ;
Sollich98 ; FieSolCat00 ; PouGut96 ; BehBiChaHenHar08 ; RedForPou11 .

###### Acknowledgements.

We thank E. Bertin for useful discussions. This work was supported by a grant
from the Simons Foundation (#454933, LB) and by a Visiting Professorship from
the Leverhulme Trust (VP1-2019-029, LB). This project has received funding
from the European Union’s Horizon 2020 research and innovation programme under
Marie Skłodowska-Curie Grant Agreement No. 893128 (RM).

## References

  * (1) Vicsek T, Czirók A, Ben-Jacob E, Cohen I, Shochet O.  Novel Type of Phase Transition in a System of Self-Driven Particles.  Phys Rev Lett. 1995 Aug;75:1226-9.  Available from: https://link.aps.org/doi/10.1103/PhysRevLett.75.1226.
  * (2) Marchetti MC, Joanny JF, Ramaswamy S, Liverpool TB, Prost J, Rao M, et al.  Hydrodynamics of soft active matter.  Rev Mod Phys. 2013 Jul;85:1143-89.  Available from: https://link.aps.org/doi/10.1103/RevModPhys.85.1143.
  * (3) Bechinger C, Di Leonardo R, Löwen H, Reichhardt C, Volpe G, Volpe G.  Active particles in complex and crowded environments.  Rev Mod Phys. 2016 Nov;88:045006.  Available from: https://link.aps.org/doi/10.1103/RevModPhys.88.045006.
  * (4) Solon AP, Fily Y, Baskaran A, Cates ME, Kafri Y, Kardar M, et al.  Pressure is not a state function for generic active fluids.  Nature Physics. 2015;11(8):673-8.
  * (5) Berthier L, Kurchan J.  Non-equilibrium glass transitions in driven and active matter.  Nature Physics. 2013;9(5):310-4.
  * (6) Berthier L.  Nonequilibrium Glassy Dynamics of Self-Propelled Hard Disks.  Phys Rev Lett. 2014 Jun;112:220602.  Available from: https://link.aps.org/doi/10.1103/PhysRevLett.112.220602.
  * (7) Ni R, Stuart MAC, Dijkstra M.  Pushing the glass transition towards random close packing using self-propelled hard spheres.  Nat Commun. 2013;4(1):1-7.
  * (8) Szamel G, Flenner E, Berthier L.  Glassy dynamics of athermal self-propelled particles: Computer simulations and a nonequilibrium microscopic theory.  Phys Rev E. 2015 Jun;91:062304.  Available from: https://link.aps.org/doi/10.1103/PhysRevE.91.062304.
  * (9) Mandal R, Bhuyan PJ, Rao M, Dasgupta C.  Active fluidization in dense glassy systems.  Soft Matter. 2016;12:6268-76.  Available from: http://dx.doi.org/10.1039/C5SM02950C.
  * (10) Klongvessa N, Ginot F, Ybert C, Cottin-Bizonne C, Leocmach M.  Active Glass: Ergodicity Breaking Dramatically Affects Response to Self-Propulsion.  Phys Rev Lett. 2019 Dec;123:248004.  Available from: https://link.aps.org/doi/10.1103/PhysRevLett.123.248004.
  * (11) Berthier L, Flenner E, Szamel G.  Glassy Dynamics in Dense Systems of Active Particles.  The Journal of Chemical Physics. 2019 May;150(20):200901.
  * (12) Janssen LMC.  Active Glasses.  Journal of Physics: Condensed Matter. 2019 Dec;31(50):503002.
  * (13) Mandal R, Bhuyan PJ, Chaudhuri P, Dasgupta C, Rao M.  Extreme Active Matter at High Densities.  Nature Communications. 2020 Dec;11(1):2581.
  * (14) Mandal R, Sollich P.  Shear-Induced Orientational Ordering in an Active Glass Former.  Proc Natl Acad Sci. 2021 Sep;118(39):e2101964118.  Available from: http://www.pnas.org/lookup/doi/10.1073/pnas.2101964118.
  * (15) Keta YE, Jack RL, Berthier L.  Disordered Collective Motion in Dense Assemblies of Persistent Particles.  Physical Review Letters. 2022 Jul;129(4):048002.
  * (16) Dunkel J, Heidenreich S, Drescher K, Wensink HH, Bär M, Goldstein RE.  Fluid Dynamics of Bacterial Turbulence.  Physical Review Letters. 2013 May;110(22):228102.
  * (17) Vicsek T, Zafeiris A.  Collective Motion.  Physics Reports. 2012 Aug;517(3-4):71-140.
  * (18) Cavagna A, Giardina I, Grigera TS.  The Physics of Flocking: Correlation as a Compass from Experiments to Theory.  Physics Reports. 2018 Jan;728:1-62.
  * (19) Szamel G.  Self-propelled particle in an external potential: Existence of an effective temperature.  Phys Rev E. 2014 Jul;90:012111.  Available from: https://link.aps.org/doi/10.1103/PhysRevE.90.012111.
  * (20) Fodor É, Nardini C, Cates ME, Tailleur J, Visco P, van Wijland F.  How Far from Equilibrium Is Active Matter?  Physical Review Letters. 2016 Jul;117(3):038103.
  * (21) Mandal R, Sollich P.  How to Study a Persistent Active Glassy System.  Journal of Physics: Condensed Matter. 2021 May;33(18):184001.
  * (22) Maloney CE, Lemaître A.  Amorphous systems in athermal, quasistatic shear.  Phys Rev E. 2006 Jul;74:016118.  Available from: https://link.aps.org/doi/10.1103/PhysRevE.74.016118.
  * (23) Flenner E, Szamel G, Berthier L.  The nonequilibrium glassy dynamics of self-propelled particles.  Soft matter. 2016;12(34):7136-49.
  * (24) Berthier L, Flenner E, Szamel G.  How active forces influence nonequilibrium glass transitions.  New Journal of Physics. 2017;19(12):125006.
  * (25) Henkes S, Kostanjevec K, Collinson JM, Sknepnek R, Bertin E.  Dense Active Matter Model of Motion Patterns in Confluent Cell Monolayers.  Nature Communications. 2020 Dec;11(1):1405.
  * (26) Szamel G, Flenner E.  Long-ranged velocity correlations in dense systems of self-propelled particles.  Europhysics Letters. 2021;133(6):60002.
  * (27) Mandal R, Sollich P.  Multiple Types of Aging in Active Glasses.  Physical Review Letters. 2020;125(21):218001.
  * (28) Liao Q, Xu N.  Criticality of the Zero-Temperature Jamming Transition Probed by Self-Propelled Particles.  Soft Matter. 2018;14(5):853-60.
  * (29) Villarroel C, Düring G.  Critical Yielding Rheology: From Externally Deformed Glasses to Active Systems.  Soft Matter. 2021;17(43):9944-9.
  * (30) Morse PK, Roy S, Agoritsas E, Stanifer E, Corwin EI, Manning ML.  A Direct Link between Active Matter and Sheared Granular Systems.  Proceedings of the National Academy of Sciences. 2021 May;118(18).
  * (31) Sethna JP, Dahmen K, Kartha S, Krumhansl JA, Roberts BW, Shore JD.  Hysteresis and hierarchies: Dynamics of disorder-driven first-order phase transformations.  Phys Rev Lett. 1993 May;70:3347-50.
  * (32) Liu C, Ferrero EE, Puosi F, Barrat JL, Martens K.  Driving Rate Dependence of Avalanche Statistics and Shapes at the Yielding Transition.  Physical Review Letters. 2016 Feb;116(6):065501.
  * (33) Ozawa M, Berthier L, Biroli G, Rosso A, Tarjus G.  Random Critical Point Separates Brittle and Ductile Yielding Transitions in Amorphous Materials.  Proceedings of the National Academy of Sciences. 2018 Jun;115(26):6656-61.
  * (34) Angelini TE, Hannezo E, Trepat X, Marquez M, Fredberg JJ, Weitz DA.  Glass-like Dynamics of Collective Cell Migration.  Proceedings of the National Academy of Sciences. 2011 Mar;108(12):4714-9.
  * (35) Kumar N, Soni H, Ramaswamy S, Sood AK.  Flocking at a distance in active granular matter.  Nature Communications. 2014;5(1):4688.
  * (36) Briand G, Dauchot O.  Crystallization of Self-Propelled Hard Discs.  Phys Rev Lett. 2016 Aug;117:098004.
  * (37) Berthier L, Biroli G, Bouchaud JP, Cipelletti L, van Saarloos W.  Dynamical heterogeneities in glasses, colloids, and granular media. vol. 150.  OUP Oxford; 2011.
  * (38) Keys AS, Hedges LO, Garrahan JP, Glotzer SC, Chandler D.  Excitations Are Localized and Relaxation Is Hierarchical in Glass-Forming Liquids.  Phys Rev X. 2011 Nov;1:021013.  Available from: https://link.aps.org/doi/10.1103/PhysRevX.1.021013.
  * (39) Berthier L, Biroli G.  Theoretical perspective on the glass transition and amorphous materials.  Rev Mod Phys. 2011 Jun;83:587-645.  Available from: https://link.aps.org/doi/10.1103/RevModPhys.83.587.
  * (40) Fily Y, Henkes S, Marchetti MC.  Freezing and Phase Separation of Self-Propelled Disks.  Soft Matter. 2014;10(13):2132-40.
  * (41) Digregorio P, Levis D, Suma A, Cugliandolo LF, Gonnella G, Pagonabarraga I.  Full Phase Diagram of Active Brownian Disks: From Melting to Motility-Induced Phase Separation.  Physical Review Letters. 2018 Aug;121(9):098003.
  * (42) ALGLIB. Unconstrained optimization: L-BFGS and CG; 2021.  Accessed: 2021-04-02.  Available from: https://www.alglib.net/optimization/lbfgsandcg.php.
  * (43) Nishikawa Y, Ozawa M, Ikeda A, Chaudhuri P, Berthier L.  Relaxation Dynamics in the Energy Landscape of Glass-Forming Liquids.  Physical Review X. 2022 Apr;12(2):021001.
  * (44) Maloney CE.  Correlations in the Elastic Response of Dense Random Packings.  Physical Review Letters. 2006 Jul;97(3):035503.
  * (45) Slaughter WS.  The Linearized Theory of Elasticity.  Boston, MA: Birkhäuser Boston; 2002.
  * (46) Mizuno H, Shiba H, Ikeda A.  Continuum Limit of the Vibrational Properties of Amorphous Solids.  Proceedings of the National Academy of Sciences. 2017 Nov;114(46).
  * (47) Lemaître A.  Structural Relaxation Is a Scale-Free Process.  Physical Review Letters. 2014 Dec;113(24):245702.
  * (48) Lerbinger M, Barbot A, Vandembroucq D, Patinet S.  Relevance of Shear Transformations in the Relaxation of Supercooled Liquids.  Phys Rev Lett. 2022 Oct;129:195501.  Available from: https://link.aps.org/doi/10.1103/PhysRevLett.129.195501.
  * (49) Karimi K, Ferrero EE, Barrat JL.  Inertia and Universality of Avalanche Statistics: The Case of Slowly Deformed Amorphous Solids.  Physical Review E. 2017 Jan;95(1):013003.
  * (50) Maloney C, Lemaître A.  Subextensive Scaling in the Athermal, Quasistatic Limit of Amorphous Matter in Plastic Shear Flow.  Physical Review Letters. 2004 Jul;93(1):016001.
  * (51) Maloney C; 2022.  Private conversation.
  * (52) Binder K, Kob W.  Glassy Materials and Disordered Solids: An Introduction to Their Statistical Mechanics.  World Scientific Publishing Company; 2005.
  * (53) Heussinger C, Berthier L, Barrat JL.  Superdiffusive, heterogeneous, and collective particle motion near the fluid-solid transition in athermal disordered materials.  EPL (Europhysics Letters). 2010;90(2):20005.
  * (54) Kob W, Donati C, Plimpton SJ, Poole PH, Glotzer SC.  Dynamical Heterogeneities in a Supercooled Lennard-Jones Liquid.  Physical Review Letters. 1997 Oct;79(15):2827-30.
  * (55) Chaudhuri P, Berthier L, Kob W.  Universal Nature of Particle Displacements close to Glass and Jamming Transitions.  Phys Rev Lett. 2007 Aug;99:060604.  Available from: https://link.aps.org/doi/10.1103/PhysRevLett.99.060604.
  * (56) Berthier L, Chandler D, Garrahan JP.  Length Scale for the Onset of Fickian Diffusion in Supercooled Liquids.  Europhysics Letters (EPL). 2005 Feb;69(3):320-6.
  * (57) Berthier L.  Time and Length Scales in Supercooled Liquids.  Physical Review E. 2004 Feb;69(2):020201.
  * (58) Henkes S, Brito C, Dauchot O.  Extracting Vibrational Modes from Fluctuations: A Pedagogical Discussion.  Soft Matter. 2012;8(22):6092.
  * (59) Bottinelli A, Silverberg JL.  How to: Using Mode Analysis to Quantify, Analyze, and Interpret the Mechanisms of High-Density Collective Motion.  Frontiers in Applied Mathematics and Statistics. 2017 Dec;3:26.
  * (60) Pal A, Kostinski S, Reuveni S.  The Inspection Paradox in Stochastic Resetting.  Journal of Physics A: Mathematical and Theoretical. 2022 Jan;55(2):021001.
  * (61) Shiba H, Kawasaki T, Onuki A.  Relationship between Bond-Breakage Correlations and Four-Point Correlations in Heterogeneous Glassy Dynamics: Configuration Changes and Vibration Modes.  Physical Review E. 2012 Oct;86(4):041504.
  * (62) Lemaître A, Caroli C.  Plastic response of a two-dimensional amorphous solid to quasistatic shear: Transverse particle diffusion and phenomenology of dissipative events.  Phys Rev E. 2007 Sep;76:036104.  Available from: https://link.aps.org/doi/10.1103/PhysRevE.76.036104.
  * (63) Guiselin B, Scalliet C, Berthier L.  Microscopic origin of excess wings in relaxation spectra of supercooled liquids.  Nature Physics. 2022;18(4):468-72.
  * (64) Lačević N, Starr FW, Schrøder TB, Glotzer SC.  Spatially Heterogeneous Dynamics Investigated via a Time-Dependent Four-Point Density Correlation Function.  The Journal of Chemical Physics. 2003 Oct;119(14):7372-87.
  * (65) Scalliet C, Guiselin B, Berthier L.  Thirty Milliseconds in the Life of a Supercooled Liquid.  Phys Rev X. 2022 Dec;12:041028.  Available from: https://link.aps.org/doi/10.1103/PhysRevX.12.041028.
  * (66) Keys AS, Hedges LO, Garrahan JP, Glotzer SC, Chandler D.  Excitations Are Localized and Relaxation Is Hierarchical in Glass-Forming Liquids.  Phys Rev X. 2011 Nov;1:021013.  Available from: https://link.aps.org/doi/10.1103/PhysRevX.1.021013.
  * (67) Caprini L, Marini Bettolo Marconi U, Puglisi A.  Spontaneous Velocity Alignment in Motility-Induced Phase Separation.  Phys Rev Lett. 2020 Feb;124:078001.  Available from: https://link.aps.org/doi/10.1103/PhysRevLett.124.078001.
  * (68) Illing B, Fritschi S, Kaiser H, Klix CL, Maret G, Keim P.  Mermin–Wagner Fluctuations in 2D Amorphous Solids.  Proceedings of the National Academy of Sciences. 2017 Feb;114(8):1856-61.
  * (69) Wysocki A, Winkler RG, Gompper G.  Cooperative Motion of Active Brownian Spheres in Three-Dimensional Dense Suspensions.  EPL (Europhysics Letters). 2014 Feb;105(4):48004.
  * (70) Nicolas A, Martens K, Barrat JL.  Rheology of Athermal Amorphous Solids: Revisiting Simplified Scenarios and the Concept of Mechanical Noise Temperature.  Europhys Lett. 2014;107:44003.
  * (71) Sollich P, Lequeux F, Hébraud P, Cates ME.  Rheology of Soft Glassy Materials.  Phys Rev Lett. 1997;78:2020-3.
  * (72) Sollich P.  Rheological Constitutive Equation for a Model of Soft Glassy Materials.  Phys Rev E. 1998;58:738-59.
  * (73) Fielding SM, Sollich P, Cates ME.  Aging and Rheology in Soft Materials.  J Rheol. 2000;44(2):323-69.
  * (74) Pouliquen O, Gutfraind R.  Stress Fluctuations and Shear Zones in Quasistatic Granular Chute Flows.  Phys Rev E. 1996 Jan;53(1):552-61.  Available from: https://link.aps.org/doi/10.1103/PhysRevE.53.552.
  * (75) Behringer RP, Bi D, Chakraborty B, Henkes S, Hartley RR.  Why Do Granular Materials Stiffen with Shear Rate? Test of Novel Stress-Based Statistics.  Phys Rev Lett. 2008 Dec;101(26):268301.  Available from: https://link.aps.org/doi/10.1103/PhysRevLett.101.268301.
  * (76) Reddy KA, Forterre Y, Pouliquen O.  Evidence of Mechanically Activated Processes in Slow Granular Flows.  Phys Rev Lett. 2011 Mar;106(10):108301.  Available from: https://link.aps.org/doi/10.1103/PhysRevLett.106.108301.