Extreme active matter at high densities

We study the remarkable behaviour of dense active matter comprising self-propelled particles at large Péclet numbers, over a range of persistence times, from τp → 0, when the active fluid undergoes a slowing down of density relaxations leading to a glass transition as the active propulsion force f reduces, to τp → ∞, when as f reduces, the fluid jams at a critical point, with stresses along force-chains. For intermediate τp, a decrease in f drives the fluid through an intermittent phase before dynamical arrest at low f. This intermittency is a consequence of periods of jamming followed by bursts of plastic yielding associated with Eshelby deformations. On the other hand, an increase in f leads to an increase in the burst frequency; the correlated plastic events result in large scale vorticity and turbulence. Dense extreme active matter brings together the physics of glass, jamming, plasticity and turbulence, in a new state of driven classical matter.

Extreme active matter, an assembly of self-propelled particles with large persistence time τp and high Péclet number, exhibits remarkable behaviour at high densities. As τp → 0, the assembly undergoes a gradual slowing down of density relaxations, as one reduces the active propulsion force f , until at the glass transition, the relaxation times diverge. In the other limit, τp → ∞, the fluid jams on lowering f , at a critical threshold f * (∞), with stresses concentrated along force-chains. As one moves away from this jamming threshold, the force-chains dynamically remodel, and the lifetime of the force-balanced configurations diverges as one approaches f * (∞), by tuning τp. In between these limits, the approach to dynamical arrest at low f , goes through a phase characterised by intermittency in the kinetic energy. This intermittency is a consequence of long periods of jamming followed by bursts of plastic yielding associated with Eshelby deformations, akin to the response of dense amorphous solids to an externally imposed shear. The frequency of these plastic bursts increases as one moves towards the intermittent phase-fluid boundary, where the correlated plastic events result in large scale vorticity and turbulence. Dense extreme active matter brings together the physics of glass, jamming, plasticity and turbulence, in a new state of driven classical matter.
Active Matter, where each particle comprising the system is driven by an internal energy source and dissipates it in movement, constitutes a new class of driven non-equilibrium systems [1]. Extreme active matter, where the magnitude of propulsion force is much higher than temperature and the direction of propulsion force persists over long times, is an extreme realisation of activity. In this limit, active systems must show strong departures from equilibrium; this expectation is borne out in active Ornstein-Uhlenbeck particles (AOUPs) at low densities [2], where steady states with finite energy-current manifest when the persistence time is sufficiently large. Even so, one might suspect that at very high densities, these distinguishing effects of activity will be firmly suppressed [3][4][5][6][7][8][9][10][11][12][13][14]. In this manuscript, we see that extreme active matter at high densities is a fount of surprises, bringing together the physics of glass, jamming, plasticity and turbulence, in a new state of driven classical matter.
Our model dense extreme active matter is a dense assembly of self-propelled soft particles subject to a propulsion force of magnitude f and whose orientation persists over a time τ p . Our main results are summarised in Fig. 5: (i) For small values of τ p , the assembly smoothly transforms from a fluid at high f to a dynamically arrested glass at low f . The phase boundary is well described by an active generalisation of RFOT theory (ARFOT) [15], with an "effective temperature" that goes as f 2 τ p /(1 + Bτ p ). However, we find that the mean kinetic energy has a different scaling behaviour with τ p , pointing to the feature that active systems should be characterised by many different "temperatures". (ii) At intermediate values of τ p , the fluid abruptly transforms into an intermittent fluid as f is lowered to f * (τ p ), characterised by intermittency in the kinetic energy. The intermittency increases as f is reduced, until at low enough f , the assembly undergoes complete dynamical arrest. (iii) This intermittency is a consequence of periods of jamming followed by bursts of plastic yielding. We identify isolated plastic events with Eshelby deformations, akin to the response of dense amorphous solids to an externally imposed shear. (iv) As one approaches f * (τ p ) from below, the plastic events become numerous and correlated in space-time, with an avalanchelike scale invariant statistics, that results in vorticity and turbulence, characterised by an inverse cascade with a Kolmogorov exponent. (v) The accumulated yielding over a time window involves the cooperative movement of a finite fraction of particles, and should manifest as a viscoelastic fluid at large time scales. (vi) In the limit τ p → ∞, the fluid reaches a jammed state on lowering f to f * (∞), with stresses concentrated along force-chains. As one moves away from this jamming threshold, by tuning τ p , the force-chains dynamically remodel, and the lifetime of the force-balanced configurations diverges as one approaches f * (∞), with an exponent z ≈ 0.71.
The stochastic dynamics of a 2-dimensional assembly of interacting active brownian particles [16], each of mass m and driven by a stochastic self-propulsion force f = f n whose direction n ≡ (cos θ, sin θ) undergoes rotational diffusion, is given by, where the i-th particle is subject to a friction γ and a thermal noise ϑ i with zero mean and variance 2k B T γδ(t − t ) that obeys fluctuation-dissipation relation. The rotational diffusion of the orientation of the propulsion force θ i is described by an athermal noise ξ i , with zero mean and correlation ξ i (t)ξ j (t ) = 2τ −1 p δ ij δ(t − t ). Its effect on the x i -dynamics is as an exponentially correlated vectorial noise with correlation time τ p , which being unrelated to the drag γ, violates fluctuation-dissipation relation. The inter-particle force f ij is modelled via Lennard-Jones interaction with particle diameter σ. See [17] for simulation details and units. Extreme dense active matter is characterised by large τ p , high Péclet number f σ/k B T and high densities. Here, we have fixed the number density to be 1.2, which is in the regime where this Kob-Andersen model [18] of passive binary soft-spheres, shows dynamical arrest at low temperatures. Further, we focus on the strictly athermal limit T = 0; we have checked that our results hold when T is small.
The set of equations (Eq. 1) for the assembly of particles are numerically integrated, using velocity Verlet algorithm, and we monitor the dynamics of density relaxations and time series of energies, stresses etc., by changing f at different values of τ p .

Low persistence time: dynamical arrest
At high propulsion force, f , the material is a fluid with time-correlations of density fluctuations measured via the self-overlap function Q(t) relaxing diffusively (see Fig. 1a(ii) and Fig. S2 in Supplementary Information [17]). As f is reduced, density relaxations slow down, until the onset of glass transition at f = f c (τ p ), estimated by fitting the variation of the relaxation times τ α versus f (Fig. 1a(iii)) with a diverging power-law [3]. The glass transition boundary, obtained for a range of τ p in this low persistence limit, can be fairly accurately described by an active generalisation of the well known Random First Order Transition (RFOT) theory [19], with an "effective temperature" that goes as Af 2 τ p /(1 + Bτ p ) (see Fig. 5), A and B being fit parameters [15]. Indeed, recent studies [3-6, 12, 13] We see that in the small t end of this log-plot, κ(t) increases linearly as t decreases and should therefore diverge, when extrapolated to t → 0. The dynamical order parameter is measured from the value of κ(t) at the earliest time that we can evaluate, i.e. t = 0 + . (Inset) Variation of the dynamical order parameter, the excess kurtosis, κex(0 + ) with f . We use the point of inflection of this curve to determine the phase boundary to the intermittent phase. (iii) The fluctuation χ4(t) = Q 2 (t) − Q(t) 2 , shows a peak at a time t for different values of f . (inset) At a fixed τp = 10 4 , the value of the peak height hp increases sharply as f approaches f * (τp) from above, then reduces again. The value of f at which hp has a maximum, for different values of τp, also marks the boundary between the liquid and the intermittent phase (see Fig. 5).
are consistent with the predictions of this active RFOT theory, in the limit of low τ p . This slowing down of particle motion, is also apparent in the time series of the mean kinetic energy E(t) as one reduces f , i.e., the mean and variance of the kinetic energy decrease as f → f c (τ p ) ( Fig. 1a(i)). However, the mean kinetic energy appears to have a nontrivial scaling with τ p (see Fig. S3 in Supplementary Information [17]), over the range of τ p values investigated. This deviation from equipartition, may be due to the fact that the joint distribution of velocities and positions at steady state do not decouple [2].

Intermediate persistence time: intermittent jamming and plastic yielding
At intermediate persistence times τ p 10 3 , the relaxation dynamics is fundamentally different from that at low τ p . At high propulsion force, f , the particles move as a fluid as before, with the time series of the mean kinetic energy showing typical Brownian fluctuations ( Fig. 1b(i)). As f is reduced, the local particle displacements start to show spatial correlations, with a growing correlation length as one approaches f = f * (τ p ) (see Fig. S4 in Supplementary  Information [17]). At and below the transition f * , the average kinetic energy (E(t)) shows sudden bursts (over a time interval τ 1 ) with periods of quiescence or jamming (over a time interval τ 2 ), typical of intermittency [20], as shown in Fig. 1b(i), characterised by large fluctuations.
To describe the dynamics of such statistical quantities that alternate between periods of quiescence and large changes over very short times, we monitor their time-dependent 4 th -moment or kurtosis [20,21]. At f > f * (τ p ), κ(t) is nearly flat and close to 3, indicating that the fluctuations are close to Gaussian. For This implies that the deformation associated with a single, isolated yield event is an Eshelby deformation.
shows an increase at small t, that becomes more pronounced with decreasing propulsion force f . We observe that κ(t) shows a power law divergence at small t, a characteristic signature of intermittency [20,21]. This allows us to describe the intermittent phase by a dynamical order parameter, the excess kurtosis κ ex (0 + ), of the increment in the kinetic energy over an infinitesimal time interval, which goes from 0 (Gaussian distribution) when f > f * (τ p ) to a finite value (indicative of broad non-Gaussian distributions) across a continuous transition at f = f * (τ p ) ( Fig. 1b(ii)). The change in this order parameter is sharp for large values of τ p (inset of Fig. 1b(ii)) and becomes more gradual as τ p is reduced, indicating that at lower τ p , the transition from liquid to intermittent phase is more like a crossover. From the variation of this dynamical order parameter over the {τ p , f } plane, we plot the non-equilibrium phase boundary in Fig. 5.
Other quantities begin to show a broad distribution as f → f * (τ p ), such as in the time-correlation of density fluctuations, Q(t). We see this in the fluctuations, χ 4 (t) = Q 2 (t) − Q(t) 2 , a measure of the dynamical heterogeneity. At fixed f , χ 4 (t) typically shows a peak at a time t (which is less than τ p ), as shown in Fig. 1b(iii), and we denote the peak height as h p . The value of h p increases sharply as f approaches f * (τ p ) from above, then reduces again (see inset of Fig. 1b(iii)). The value of f at which h p has a maximum, for different values of τ p (see Fig. S5 in Supplementary  Information [17]), provides another marker of the boundary between the liquid and the intermittent phase, indicated by blue circles in Fig. 5.
Within the intermittent phase, we notice that the sudden increase in kinetic energy during a burst is instantaneously accompanied by a non-reversible release in the potential energy ( Fig. 2a(ii), as well as visible spikes in the local shear stress (Fig. 2a(iii)). Thus, driven by persistent active stresses, configurations of particles in the intermittent phase experience a buildup of the elastic stress, a transient jamming (E = 0), followed by sudden yielding, seen as a burst of kinetic energy ( Supplementary Movies [17]).
The bursts in kinetic energy are accompanied by local structural reorganisations associated with sudden collective non-affine displacements of a finite fraction of particles (see Supplementary Movie 1 [17]). These bursty features are apparent in the thresholded displacements over a time τ (see Fig. S6 in Supplementary Information [17]), and more directly in the spatial maps of the particle displacements. This implies that the intermittent steady state exhibits a continual yielding and jamming of macroscopically large structures.

Plastic yielding: Eshelby deformations
Deep in the intermittent phase, these bursts in kinetic energy and associated plastic yielding, are rare and isolated (Fig. 2a), allowing us to analyse the deformation field around a single burst event.
The radial and azimuthal components of the displacement field u(r, θ) surrounding the single yielding event, shows a quadripolar symmetry (Fig. 2b(i)-(iii)) and a long range decay with radial distance that goes as 1/r (Fig. 2b(iv)). A similar feature is shown by the local elastic shear stress propagated as a consequence of a single yielding event. This is the well studied Eshelby deformation profile [22], which describes elementary local deformations in an amorphous solid under external uniform shear [23,24]. The unexpected appearance here of the Eshelby stress, is a result of local shear arising from internal stirring at the scale of the active particle.
As one moves towards the intermittent phase -liquid boundary from below, the bursts get more frequent and are bunched up. The distribution of the periods of intermittent bursts (τ 1 ) and quiescence (τ 2 ), is a power law with an exponential cut-off: P (τ 1 ) ∼ τ −α 1 exp(−τ 1 /τ 10 ) and P (τ 2 ) ∼ τ −β 2 exp(−τ 2 /τ 20 ) (see Fig. 3a for fit parameters), with the cut-off moving to larger times as f → f * (τ p ) from below. In the vicinity of f * (τ p ), the distributions are power-laws, P (τ 1 ) ∼ τ −2.6 1 and P (τ 2 ) ∼ τ −1.85 2 (Fig. 3b). Each of these plastic events give rise to stresses that propagate through the material (see Supplementary  Movies 2&3 [17]). Outside the plastic zones, the rest of the material should respond elastically, with an anisotropic elastic kernel. However, the occurrence of multiple yielding events will result in strong correlations between events, one triggering another, that will make the kernel isotropic, since the directions of local shear due to active forcing would be randomly oriented.

Accumulated yielding, turbulence
Since, for a small enough f , the plastic bursts are bunched up discrete events, the number of particles n c (∆t) that undergo irreversible displacement within a time window ∆t is the number of displaced particles per event times the number of events within the window ∆t. We find that n c ∼ N 2 (Fig. 3c), suggesting that each intermittent yielding event, involves the collective displacement of a finite fraction of particles [25]. This implies that the occurrence of more and more such events will cause the material to flow at long time scales, with a time dependent viscosity η(t), determined by the accumulated yielding upto time t, eventually reaching a constant steady-state value for t τ p . This shows up, at the level of single tagged particle dynamics, as eventual diffusive motion when t τ p (see Fig. S7 and Fig. S9 in Supplementary Information [17]), and coincides with the relaxation of the self-overlap function Q(t) [17].
As we reduce the active force f , the tagged-particle diffusion coefficient decreases (Fig. S7 in Supplementary Information [17]), and eventually vanishes as one approaches dynamical arrest. As before, we obtain the dynamical arrest boundary f = f c (τ p ) (see Fig. 5), by fitting the data for the α−relaxation time τ α measured from Q(t) (Fig. S8 in Supplementary Information [17]), to a power-law divergence.
On the other hand, as we increase the forcing f towards the phase boundary f * (τ p ), we observe that the scale-free intermittency displays a kind of plastic turbulence [26,27] in an actively stirred dense material. This is seen in the spatiotemporal dynamics of the displacement fields which show large swirls; see Supplementary Movie 4 [17], and in the power-spectrum of the kinetic energy density. The intermittent jamming-yielding due to local active stirring transfers energy from small scales and dissipates it over larger scales, leading to an inverse cascade, where the energy spectrum crosses over from E(k) ∼ k −5 to E(k) ∼ k −5/3 at lower k (Fig. 3d). The crossover to the Kolmogorov spectrum happens at a scale corresponding to the scale of vorticity [28]. This stress production and dissipation gives rise to a non-equilibrium steady state with a finite energy-current.

Infinite persistence: jamming/unjamming
Analysis of the τ p = ∞ limit, brings in a new facet of extreme active matter. This limit corresponds to a situation where the initial directions of particle self-propulsion are quenched in random directions. From being a fluid with mobile particles at large f , the assembly jams at f * (∞) = 1.67, where the kinetic energy goes to zero as ∼ |f − f * (∞)| 3/2 (Fig. 4a,b). The distribution of forces P (F) changes from a broad distribution with exponential tails to a deltafunction at F = 0 at f * (∞); the jamming transition is associated with a force-balanced configuration of the soft particles (Fig. 4c). The width of P (F) and the mean kinetic energy go continuously to zero as f approaches f * (∞) (inset Fig. 4c). This allows us to identify f * (∞) as a jamming/unjamming force threshold for active yielding. As discussed in a recent work [29], the density-dependent f * (ρ, ∞) will trace out a yielding line in the jamming phase diagram [30] of dense amorphous materials, with active forcing being the control variable. Thus, we expect to find critical behaviour in the proximity of f * (∞).
In the vicinity of f * (∞), we map the contact network of particles and evaluate the total force between pairs of soft particles in contact; we find that these forces are distributed along force-chains (Fig. 4d(iii)). With increasing f , the force chains dynamically remodel, as the structures relax; see Fig. 4d(i),(ii). Likewise, dynamical reorganisation of the force network also occurs when we move slightly away from the jammed regime, by decreasing τ p , keeping f . The data is shown for f = 1.6, τp = 10 4 (marked by red triangle in Fig. 5). (c) The intermittent yielding events involve the non-affine displacement of a finite fraction of particles, as seen in this plot of nc, the number of particles that show a non-affine displacement within a time window ∆t = 10 4 , versus total particle number N . (d) The scale-free intermittency close to the phase boundary f * (τp), is associated with plastic turbulence as seen in the power-spectrum of the energy density E(k) which shows an inverse cascade from an injection scale, shown by the arrow. Data shown for f = 1.6, τp = 10 4 . The crossover from a steep spectrum k −5 to the Kolmogorov spectrum k −5/3 at lower k, is set by the scale of the vorticity.
fixed. Under these conditions, the dynamics of force-chains typically show periods of jamming in a force-balanced configuration, interspersed with bursts of remodelling. As one approaches the jamming threshold at τ p = ∞, the mean lifetime of the force-balanced configurations diverges as τ F ∼ τ z p , with a new dynamical critical exponent z ≈ 0.71 (Fig. 4e).

Discussion
Putting all this together our study suggests a rather rich phase diagram (Fig. 5). The study of extreme dense active matter allows us to explore the crossover between glass physics, where the dynamics proceeds by density relaxation, and jamming-yielding physics, where the dynamics is controlled by stress buildup and release via macroscopic flows. We emphasize that the intermittent plastic deformation and turbulent flows are constitutive and not response to an externally imposed stress.
While our present study was done at T = 0, we have checked that including a small temperature via a thermal noise ϑ i gives similar results, as long as the Péclet number is high; the crossover behaviour from these different regimes are likely to be quite subtle.
We have fixed the overall density of particles in our study, however [29] shows that the athermal jamming transition at τ p = ∞ occurs over a range of densities, making this active jamming critical point, density dependent. If however we make the density very low, while still keeping τ p = ∞, we would arrive at a jammed gas phase, with isolated islands of jammed material [31].
Are there natural or synthetic realisations of extreme active matter? Herds of animals, such as penguins or bulls, dense collection of vehicles, ants or microbots and even trite examples such as a herd of rugby players, could be possible realisations. Promising candidates for extreme active matter are monolayers of persistently motile cells; in- deed [32] observe jamming-yielding behaviour in such monolayers of epithelial cells. It would be challenging however to construct synthetic realisations of extreme active matter, and we eagerly look forward to controlled experimental studies on these. Dynamical phase diagram in f − τp at fixed density. At low τp (0 < τp < 10), there is a direct transition at f = fc (red dots), from the liquid to a dynamically arrested phase, as measured from the the divergence of τα, the slow relaxation of the density fluctuations ( Fig. 1a(i)). This transition to the dynamically arrested state continues into larger values of τp (red dots). The dashed line is a fit to the active RFOT theory [15] with an effective temperature Af 2 τp/(1 + Bτp), where A and B are fit parameters. At larger values of τp > 10, a new intermittent phase, intervenes between the liquid and the dynamically arrested phase. The liquid-intermittent phase boundary is obtained from the dynamical order parameter, the excess kurtosis κex (Fig. 1b(ii)) and the peak height hp (Fig. 1b(iii)), shown with blue circles. This transition is very sharp and continuous at large values of τp, and gets progressively broader and less defined at lower values of τp (colour bar represents the value of κex(0 + ) in log 10 for different f and τp). This intermittent phase is characterised by bursts of plastic yielding and plastic turbulence close to the transition to the liquid phase. At τp = ∞, the assembly shows a sudden transition from a liquid to a jammed configuration at a force threshold, where the particles suddenly get into a force-balanced configuration (Fig. 4d).

A. Model and Methods
For our study, we consider a well-studied two-dimensional model, viz., a 65:35 binary Lennard-Jones (LJ) mixture [18], with particles interacting via the following potential, where r is the distance between the i-th and the j-th particle, i.e., r = |r i − r j | and α, β represent either A-type or B-type particles. The strength and the range of the interaction are set by αβ and σ αβ respectively. In our simulation we have chosen the values of σ αβ and αβ to be: σ AB = 0.8σ AA , σ BB = 0.88σ AA , AB = 1.5 AA , BB = 0.5 AA . The composition of the A : B mixture helps to avoid crystallisation, in the absence of activity. The potential has been truncated at r c αβ = 2.5σ αβ and has been shifted accordingly such that both the potential and the force remain continuous at the cut-off. The unit of length and energy in our simulation are set by σ AA = 1 and AA = 1 and the study is done for an overall number density of ρ = 1.2. All the particles have the same mass (m = 1) and the time unit is τ LJ ≡ mσ 2 AA / AA = 1. We study the assembly under athermal conditions, with a self propulsion force of magnitude f along a unit vector associated with each particle, in the presence of a Langevin bath. Thus, the equation of motion for each particle is, where m is the mass of a particle, v i is the velocity of the i-th particle, γ is the friction coefficient. For our study, we choose γ = 1. Thus, the particle inertia relaxation time τ γ ≡ m γ is comparable to τ LJ . f ij is the LJ interaction force between the i and j-th particle, n i is the unit vector associated with the i-th particle along which the propulsion force is being imparted, f is the strength of the propulsion force. See Fig. S1 for a schematic snapshot of a typical configuration, with the propulsion direction of each particle being indicated by an arrow. The dynamics of the direction of n i follows simple rotational diffusion equation with diffusion constant D R ∝ τ −1 p (τ p is the persistence time). In our study, we will tune both f , and τ p , and study the dynamical behaviour of the assembly of the active particles. The number of particle used in the simulation varies between N = 1000 − 10000. All data presented here have also been averaged over 32 − 96 independent realisations, if not mentioned otherwise. To characterise the dynamics of the dense liquid, under steady state or transient conditions, as we vary the active forcing (f ), for various choice of τ p , we measure the mean squared displacement (M SD) ∆ 2 (t) and self-part of the two-point overlap correlation function Q(t), defined as, where, · · · represents an average over the time origin t 0 , N is the number of particles in the system and the parameter b is associated with the typical vibrational amplitude of the caged particles. Throughout our analysis we have used b = 0.3 and we have verified that our results are insensitive to moderate changes in b.
Plots showing the variation of ∆ 2 (t) and Q(t), with f and τ p are shown in Fig. S2, S7, S9. Another quantity of interest is the displacement field, calculated over a time scale of τ α ; see Fig. S4. Using that, we calculate the two point spatial correlation function, defined as C(r) = x τα ( r)· x τα ( 0) where x τα ( r) is the displacement of the particles at position r over the timescale τ α . To extract the lengthscale (ζ) associated with the size of this cooperatively rearranging regions we use the relation C(ζ) = 1/e, and we monitor how ζ varies with f ; see Fig. S4.    (Bottom) Corresponding mean squared displacement, ∆ 2 (t). For τp ≥ 10 3 , the characteristic relaxation timescale is t/τp ≈ 1, and diffusive motion is also seen to set in, beyond this timescale.