Active Jamming at Criticality

Jamming is ubiquitous in disordered systems, but the critical behavior of jammed solids subjected to active forces or thermal fluctuations remains elusive. In particular, while passive athermal jamming remains mean-field-like in two and three dimensions, diverse active matter systems exhibit anomalous scaling behavior in all physical dimensions. It is therefore natural to ask whether activity leads to anomalous scaling in jammed systems. Here, we use numerical and analytical methods to study systems of active, soft, frictionless spheres in two dimensions, and elucidate the universal scaling behavior that relates the excess coordination, active forces or temperature, and pressure close to the athermal jammed point. We show that active forces and thermal effects around the critical jammed state can again be captured by a mean-field picture, thus highlighting the distinct and crucial role of amorphous structure in active matter systems.

Disordered systems of particles are ubiquitous in nature, from molecular glasses, colloidal suspensions and foams to biological tissues and granular materials.As such, models of jammed solids have naturally attracted intense attention from both theorists and experimentalists alike .For an athermal system of soft frictionless spheres, it is by now well established that upon increasing its density, the system generically develops a yield stress under which the solid respond elastically [4][5][6][7][8][9][10][11][12].This jamming transition, known as "point J", is characterized by a critical packing fraction ϕ c , at which the jammed solids is marginally stable [1,4,6] and displays an isostatic contact network [6].
Intriguingly, the mechanical properties of such a system exhibit critical scaling in the vicinity of point J [27][28][29], and the corresponding upper critical dimension, d u , is believed to be 2 [4,27,[30][31][32][33][34] Importantly, the fact that d u = 2 has enabled researchers to successfully use mean-field theory to elucidate certain critical behavior of jamming for physically relevant dimensions (d = 2, 3).
Interestingly, a recent study revealed that systems of active Brownian particles (ABP) at high density jam intermittently and that the lifetimes of these transiently jammed states lengthens with the persistence of the particles [95].Since the lifetimes of these transient jammed states can be arbitrarily long, they are clearly experimentally relevant.However, the associated scaling behavior remains unexplored.In this Letter, we fill this void by elucidating the scaling behavior of these transiently jammed states under both activity and thermal fluctuations around the static critical jamming point.
Scaling ansatz.-Startingfrom the established fact that passive and athermal jamming is a critical phenomenon for which the critical point is exactly at isostaticity (i.e., the average coordination z is z iso = 2d), analysis based on scaling ansatz have provided scaling relations among diverse physical observables [27].Here, we also start by formulating a scaling ansatz that incorporates the effects of active forces in the athermal jamming scenario.Since prior work has demonstrated that active forces have the generic effects of unjamming a system, we start by considering a generic system of soft frictionless spherical particles above jamming onset, i.e. with pressure p > 0. In the absence of active forces, the following scaling relation between ∆z ≡ z − z iso and p is well established in 2D systems [27]: ∆z ∼ p 1/1.88 .
(1) We now imagine that the particles can exert a persistent active force f in a randomly chosen direction.As adding this force will tend to unjam the system, we anticipate that ∆z should decrease with f .At the same time, we expect the system to revert back to the same inherent structure when active forces are removed as long as ∆z remains positive.In other words, we expect the system to only explore locally its potential energy landscape and to remain in the basin of attraction of its initial static configuration.This expectation is justified by previous numerical simulations [95] and our own simulation results.
Around this critical region where ∆z, p, and f are all small, we expect that the system to be scale invariant with a single length scale that is controlled by one of the control parameter [96]; we choose our control parameter to be p in this case since it is the physical quantity that we actually control in most simulation studies.As such, we arrive at the following scaling ansatz: where S is a universal scaling function such that lim x→0 S(x) is a positive constant so that we recover the static scaling relation between ∆z and p (1).Again, since increasing f has the generic effect of unjamming the system, we anticipate that the scaling function S(x) is monotonically decreasing with respect to x, and eventually becomes 0 at x c .In other words, there exists a critical force f c such that when ∆z is exact zero in the system.We now test this scaling ansatz via means of numerical simulations.
Model and simulations.-Westudy the dynamics of two-dimensional (d = 2) jammed packings of N spherical particles interacting via purely repulsive harmonic forces when subjected to active forces and thermal fluctuations.We start by creating static jammed packings at fixed pressures p providing a way to control the packings distance to the athermal jamming transition point.To do so, we start with random configurations of discs at high volume fraction (ϕ = 0.95) in a square box of size L with periodic boundary conditions.We then follow the FIRE (Fast Inertial Relaxation Engine) energy minimization procedure until the maximum unbalanced force on any particle is less than f t = 10 −14 .Following each energy minimization step, we either increase (if the pressure is below target) or decrease (if the pressure if above target) the diameter of the particles proceeding to a bisection until we obtain a mechanically stable configuration within 1% of the target pressure p.To avoid crystallization in the system, we work with 50 : 50 binary mixtures of particles with diameter ratio 1 : 1.4 as is typically done in studies of disordered solids [6,8].
These mechanically stable configurations form the initial conditions of our numerical simulations.Here, the time evolution of the positions of the particles is generically governed by an overdamped Langevin equation of the form where r ij = |r i − r j | is the distance between particles i and j and êi = (sin θ i , cos θ i ) denotes the direction of self-propulsion, f is the strength of active force and γ is a friction coefficient.Both η i and ξ i are zero-mean, unit variance Gaussian random variables such that The direction of the self-propulsion is thus governed by rotational diffusion with coefficient D r .Finally, U (r ij ) denotes the pairwise particle interactions which we model through a purely repulsive linear spring potential such that where σ ij = (σ i + σ j )/2 is the average diameter of particles i and j and Θ(x) is the Heaviside function.In what follows, we compare two main scenarios: (i) subjecting static jammed packings to persistent active forces via an ABP model, corresponding to solving Eq. ( 4) with D = 0 and f > 0 and (ii) subjecting static jammed packings to thermal fluctuations via Brownian dynamics, corresponding to solving Eq. ( 4) with D > 0 and f = 0. Finally, we also perform simulations for active Ornstein Uhlenbeck particles (AOUP) particles (see [97] for details).
We nondimensionalized the equations of motion using the smaller particle diameter σ and potential energy prefactor ε as basic units of length and energy.In ABP simulations, the typical timescale is given by the persistence time τ p = D −1 r , while in Brownian dynamics simulations, a more natural choice for the timescale is given by τ = σ 2 /D.In all simulations, we set σ = 1, ε = 1 and γ = 1.Further, the rotational diffusion coefficient D r is taken to be 10 −2 unless stated otherwise.To discard any potential finite size effect in our results, we varied the number of particles in the system from N = 128 to N = 4096.The target pressures of the initial configurations were taken in the range [10 −7 , 10 −1 ].Finally, results are averages over 100 independent realizations.
To estimate the scaling exponent α using Eq. ( 3), we start with static packings above jamming onset with different pressure values p in the absence of active forces; then increasing the active force strength, we find the critical value of the active force for which the time-averaged excess coordination ∆z i becomes zero.Fig. 1 shows that α = 1.5 ± 0.1 .
We now present a heuristic argument that may justify the value of α obtained from our simulation.Heuristic argument for estimating α.-Under active forcing, the unjamming transition is generically controlled by the excess coordination ∆z; it thus makes sense to expect the value of the critical active force to be dictated by the statistics of interparticle contacts in the overjammed static configurations.The scaling behavior of the distribution of interparticle gaps h has been studied both analytically and numerically [31,34,[98][99][100][101]. Specifically, the probability distribution of h in the vicinity of the athermal jamming point is given by In our soft sphere simulations, we obtain γ = 0.42 ± 0.04 [97] which is close to the value of γ ∞ = 0.41269... predicted in the mean-field theory of the glass transition in hard spheres in the limit of infinite dimensions [31,91,100].In our model of soft spheres, the gap distribution shows a weak dependence on the pressure, which we neglect here as it does not affect our argument (at least within the precision of our numerical estimate (7)).
We now assume that the gap distribution remains the same at small enough pressures p and active forces f .As the pressure p increases in the absence of active forces, some of these gaps close and new contacts form (see cartoon in fig.2).Such an increase in new contacts is given by ∆z ∼ δ 0 h −γ dh ∼ δ 0.58 (9) where δ is the average overlap.Note that since δ is known to scale like p in the overjammed regime [100], the numerical value is already in good agreement with the known scaling law between ∆z and p (1).We now consider that each particle is subject to an active force f .To eliminate all the newly formed contact and thus bring the system back to the isostaticity, the average force needed is where the upper-bound δ again corresponds to the limit at which N ∆z contacts are eliminated.This heuristic argument thus support the value of α (7) obtained from our simulation.
Universal scaling function S.-This newly found scaling allows us to elucidate the form of the universal scaling function S numerically.As presented in Fig. 3(a), S(x) is monotonically decreasing with f for all pressures.We obtain S(x) by collapsing the surface plot displayed in [97].Provided the scaling behaviors given in Eqs. ( 1) and (3), we attempt to collapse our data.In particular, as the excess coordination scales as p ∼ ∆z 1.88 , we can use this information to scale the curves in the low activity regime.Further, using the scaling of the critical active force with pressure, we obtain a good collapse of ∆z as a function of f for all pressure values as can be seen on Fig. 3(a).This collapse of our data supports the existence of scaling relations and existence of universal behavior in the vicinity of the jamming point, i.e. at low enough pressures.
Left continuity of S as it vanishes.-Beforegoing any further, it is important to discuss one of the conditions under which measuring the critical active force f c is meaningful.In the use of the scaling relation between f c and p in Eq. ( 3) to obtain the scaling exponent α above, we have focused on the vanishing point S (or, equivalently, the point where ∆z = 0).However, it is clear that the system can no longer be rigid when ∆z becomes negative (i.e.below isostaticity) and so the scaling ansatz cannot be correct in this negative ∆z regime.In other words, it is unclear whether the location when S vanishes is continuous at all.Fortunately, left continuity of S-which corresponds to dS(x)/dx| x=xc,− being welldefined-is all that is required for the above argument to work.In particular, this condition implies the following scaling relation for where A ≡ −dS(x)/dx| x=xc,− is a positive constant.We show in Fig. 3(b) that the above scaling relation is indeed satisfied; this, in turn, implies the left-continuity of S at the vanishing point, justifying our procedure to determine α (see Fig. 1).
Scaling behavior and persistence.-Studies of athermal jammed packings generically characterize jammed configuration as having a positive excess coordination.When analyzing systems subjected to active forcing, we chose to extend this criterion and define jammed systems as systems whose steady-state time-averaged excess coordination is positive.However, nothing prevents the systems coordination from falling below the isostatic limit at least transiently.As a consequence, our systems may transiently unjam before rejamming.To analyze this, we studied (i) distributions of Debye-Waller factors, (ii) mean square displacements and (iii) intermediate scattering functions in our packings under active forces.We note that strikingly both quantities display signatures of the critical transition at f = f c and show that structural rearrangements are rare even on timescales which are long compared to the active force persistence time Having focussed on the long persistence time regime so far, we now consider how robust the scaling regime is as the persistence time of the active forcing decreases.Fig. 4(a) shows that the scaling behavior indeed persists, under the same measurement protocol, even as persistence length goes down by two orders of magnitude.Interestingly, we observe that as the persistence length of the particles decreases, the regime of validity of the critical force scaling shifts towards larger pressure values (as evidenced by the shoulders developing at the high p limit), while scaling (3) remains intact in the vicinity of the jamming point.
Active forcing vs. thermal fluctuations.-Aswe have seen that the scaling regime remains robust even for very small persistence lengths, it is thus natural to extend our analysis to the situation where jammed packings are subjected to thermal fluctuations (corresponding to the zero persistence time limit) and ask the question of whether the same scaling can also be measured robustly in this case.Although the jammed system will inevitably melt under thermal fluctuations in the long time limit, we will test the robustness of the scaling regime in the intermediate time regime (as defined within the same framework of our standard simulation and measurement protocols used so far).
To study thermal fluctuations, we switch off the activity and follow dynamically as above the network of contacts in Brownian dynamics simulations.We identify the critical temperature T c at which the time-averaged excess coordination ∆z vanishes.To be able to compare directly the effects of thermal fluctuations to that of active forces, we define here our critical forcing as f c,2 = √ 2T c .Fig. 4(b) shows the critical thermal forcing f c,2 as a function of the pressure p in systems with N = 2048 particles (see [97] for a system size analysis).Strikingly, we show that the critical forcing for both active Brownian particles and their passive counterpart show the same power-law scaling with pressure.Finally, we also verify that we obtain the same scaling for AOUPs for which we estimate the critical force as f c,3 = 2D a /τ a , where τ a is the persistence time of the active force magnitude (for details see [97]).We thus conclude to the existence of dynamically jammed states even at finite forcing whose scaling behavior is independent of the driving mechanism in the vicinity of point J.
Summary & outlook.-Insummary, we have studied the effects of active forces and thermal fluctuations on long-lived jammed states of soft frictionless particles in two dimensions close to the critical athermal jamming point.Using numerical simulations supplemented by a heuristic analytical argument, we elucidated the universal scaling relations between the excess coordination, the active force (or temperature), and pressure.In particular, the heuristic, mean-field based argument that we used to support the scaling exponent suggests that the mean-field picture applicable to athermal jamming continues to hold.This is in surprising contrast to diverse active matter systems in which anomalous scaling behavior is the norm.An interesting future direction would be to seek the existence of diverging length or time scales.
This work was supported by Leverhulme Trust (Grant: RPG-2019-374).Authors acknowledge the support of Research Computing Services at Imperial College London.static packing.Being less connected than other particles, floaters have peculiar dynamical properties.In particular, when the systems is subject to thermal or active fluctuations, they move much more than the other particles and so their dynamics, while not representative of the overall dynamics of the system, may dominate the sum in Eq. (S7).One way to deal with floaters is to remove them from the sum in Eq. (S7); to do so, we need to identify the floaters which is not an easy task in a dynamic packing.In an overcompressed static packing, floaters are identified based on force balance considerations.Namely, they are defined as particles having strictly less than 3 contacts and easily removed recursively.However, this method can not be used for either unjammed states or dynamical states (as particles may come in and out of contact).Instead, our definition of floaters follows from [2] where floaters are defined as particles whose DW factor f DW is 5 times larger than the median DW factor fDW .Once identified, the N f floaters are removed from the sum in Eq. (S7) which thus runs to N ′ = N − N f .The distributions of observed DW factors P (f DW ) in configurations subject to active forces is very informative.As seen in Fig. S5, static packings at pressure p = 10 −4 subject to active forces below the critical active forces display a bimodal DW factor distributions [see Fig. S5(a)-(c)].In these bimodal distributions, the peak at higher DW factors is much smaller in amplitude and is a clear signature of the floaters.As the driving active force is increased above the critical active force measured via the excess coordination, we observe a merging of the two peaks in the distribution which becomes unimodal [see Fig. S5(d)-(e)].This is clear in Fig. S5(f) where we represent the distribution of rescaled DW factors fDW = f DW / fDW , where fDW denotes the median of the distribution.
Armed with a definition for the floaters, we can now compute unambiguously the mean-squared displacement ∆ 2 (t) for systems subjected to active forces below and above the critical active force.Fig. S6(a) shows that for all active forces f < f c , the MSD clearly plateaus and our packings are fully caged (over the timescale of our long-time simulations).For f > f c , we observe a strikingly different behavior with diffusive MSDs and we conclude that the system is unjammed passed the critical active force f c .Finally, we compute the self intermediate scattering function (ISF) over the same simulations.The ISF is defined as where the wave vector k was chosen to be proportional to the inverse of the position of the first peak in the radial pair correlation function g(r).Fig. S6 shows clearly that the ISF does not significantly decay over the full timescale of our simulations when f < f c ; we thus conclude to the lack of structural rearrangements over the course of these simulations.On the other hand, for f > f c , we observe a very fast and sharp decay of the ISF to zero clearly denoting unjammed states.Altogether, these observations support our choice of excess coordination as a key parameter to identify the critical strength of active force for the transition.

FIG. 2 .
FIG.2.Schematics illustrating our heuristic argument for value of α.(a) Particle configuration at critical jamming under no active force.The gap distribution between nontouching neighbors is expected to follow the scaling law in(8).(b) Upon further compression (or equivalently, upon particles' expansion), the pressure p becomes positive and new contacts are formed as gaps close because the particles move closer to each other.The new particles' profiles are shown in solid lines whereas the pre-compressed profiles are in dashed lines.(c) When active forces are switched on (f > 0), contacts can be broken again as the active force (blue arrow) of the red particle can counteract the steric interactions with its neighbors.

1 FIG. 3 .
FIG. 3. Universal scaling function S(x).(a) Collapsed data for excess number of contacts per particle (∆z) as a function of the active force f at various pressure values for N = 2048.A good collapse is obtained via a rescaling of ∆z by p β (with β = 1.88 for d = 2, see [27]) and f by p α with α = 3/2.(b) Evolution of the excess coordination number ∆z near the critical active force for various values of pressure for N = 2048.The black solid line shows a linear scaling.

1 3 / 2 FIG. 4 .
FIG. 4. Effect of persistence and thermal fluctuations.(a) Critical active force as a function of pressure for various values of the rotational diffusion constant Dr in simulations of N = 512 active Brownian particles (ABP).The solid black line shows a power law scaling with exponent α = 3/2.For lower values of Dr, the self-propulsive forces become more persistent and the regime of validity of our power-law scaling is limited to lower values of the pressure.(b) Comparison of critical forcing for three perturbation protocols for systems of size N = 2048: (i) persistent active forces with constant magnitude but diffusive direction in active Brownian particles simulations (ABP), (ii) thermal fluctuations in Brownian dynamics simulations (BD) and (iii) persistent active forces with fluctuating magnitude and direction in active Ornstein-Uhlenbeck particles simulations (AOUP).The critical forcing follows in all three cases the same scaling behavior as a function of p close to the passive athermal jamming point (p → 0).The solid black line shows a power law scaling with exponent α = 3/2.
FIG. S5.Distribution of Debye-Waller (DW) factor for systems initially at pressure p = 10 −4 for strength of active forces both below and above the critical active force fc.Active force magnitudes are given as follows (a) f = 10 −9 , (b) f = 10 −8 , (c) f = 10 −7 , (d) f = 2 × 10 −6 and (e) f = 3 × 10 −6 .At pressure p = 10 −4 , the critical active force is measured to be fc = 2.24 × 10 −6 .(f) Distribution of rescaled and shifted DW factors fDW.Blue lines denote distributions with f < fc exhibiting a clearly bimodal distribution characteristic of a mechanically stable vibrating backbone accompanied by floaters while orange lines denote distributions for f > fc for which bimodality has disappeared.

F 1 FIG
FIG. S6.(a) Mean-squared displacement ∆ 2 (t) and (b) self intermediate scattering function Fs(k, t) measured in systems initially at pressure p = 10 −2 and subject to active forces below and above the critical active force fc = 0.0245.As in Fig. S5, color encodes the value of the active forces relative to the critical active force.