Dynamic heterogeneity in active glass-forming liquids is qualitatively different compared to its equilibrium behaviour

Activity driven glassy dynamics is ubiquitous in collective cell migration, intracellular transport, dynamics in bacterial and ants colonies as well as artiﬁcially driven synthetic systems such as vibrated granular materials, etc. Active glasses are hitherto assumed to be qualitatively similar to their equilibrium counterparts at a suitably deﬁned eﬀective temperature, T eff . Combining large-scale simulations with analytical mode-coupling theory for such systems, we show that, in fact, an active glass is qualitatively diﬀerent from an equilibrium glassy system. Although the relaxation dynamics can be similar to an equilibrium system at a T eff , eﬀects of activity on the dynamic heterogeneity (DH), which has emerged as a cornerstone of glassy dynamics, is quite nontrivial and complex. In particular, active glasses show dramatic growth of DH, and systems with similar relaxation time and T eff can have widely varying DH. Comparison of our non-equilibrium extended mode-coupling theory for such systems with simulation results show that the theory captures the basic characteristics of such systems. Our study raises fundamental questions on the supposedly central role of DH in controlling the relaxation dynamics in a glassy system and can have important implications even for the equilibrium glassy dynamics.

Active glass, characterizing the slow dynamics in a dense-active system of self-propelled particles (SPP), with a self-propulsion force, f 0 , and a persistence time, τ p , of their motion, remains an intriguing problem of statistical physics.A number of recent experiments show that glass-like dynamics, that is the extreme dynamical slow-down of a system without any discernible structural change or phase transition [1,2], is important for the dynamics in several biological and biology-inspired systems, such as a cellular monolayer [2,3], cell-cytoplasm [4][5][6][7], colonies of bacteria [8] and ants [9,10], vertically vibrated rods [11][12][13], and several other artificial active systems [14][15][16].These systems can be conveniently modeled via a dense active matter of SPPs [17,18].Motivated by the experimental findings, such systems have been studied via simulations [19][20][21][22][23][24] as well as extended analytical theories of equilibrium glasses, such as mode-coupling theory (MCT) [25][26][27][28][29][30] and random first order transition theory [29].Studies on single-particle dynamics [31] and jamming transition [19,32,33] have also been extended for active systems.The current consensus in the field is that activity pushes the glass and jamming transitions to lower temperature or higher density, though the nature of an active glass seems to be qualitatively similar to an equilibrium glass and the difference lies in the quantitative description, such as in the form of long-ranged velocity correlation [30,34] or the evolving effective temperature, T ef f , much like a sheared system [25,26].However, this apparent similarity of glassy dynamics in an active and equilibrium system is somewhat puzzling as activity is known to exhibit nontrivial behavior, such as flocking [35] or giant number fluctuation [17,18,36,37] in a dilute system, and calls for rigorous theoretical exploration of the problem.
Apart from developing a coherent framework for the slow dynamics in these biological systems, such theoretical studies can also provide deeper insights into the physics of equilibrium glasses as they extend the scope and extent of the problem in the form of control parameters as well as emergent behaviors.In this work, through a combination of detailed large scale molecular dynamics simulations and an analytical study via extended active MCT [26], we show that, in fact, an active glass is qualitatively different from its equilibrium counterparts: whereas the relaxation dynamics is similar to that of equilibrium glasses, dramatic effects of activity is observed in the dynamic heterogeneity (DH).DH refers to the coexistence of heterogeneous transient local dynamic behavior and has emerged as a salient feature of glassy dynamics [1,38,39].Confirmed both in simulations [40][41][42] and experiments [12,43], existence of DH in glassy systems provide the much sought after growing length scale to interpret the problem as a critical phenomena with an underlying transition where the DH length scale controls the relaxation dynamics in a glassy system.Compared to equilibrium systems, active glasses show dramatic growth in DH with increasing activity and systems with similar relaxation times can have varying DH suggesting a possible decoupling from relaxation dynamics.
Although most works till date have focussed on the relaxation dynamics in an active glass, there do exist a handfull of studies [16,[44][45][46] investigating different aspects of DH revealing nontrivial effects of activity: Avila et al found [44] strong DH and large dynamic length scale in a driven granular fluid compared to those in equilibrium systems, experiments on active cellular monolayers [45,46] show higher DH with increasing activity though the relaxation time becomes smaller, a trend opposite to what one expects for equilibrium glasses.Our theoretical study rationalizes these puzzling results.We first show that the relaxation dynamics of an active glass can be described similar to that of an equilibrium system at an effective temperature, T ef f (Fig. 1), however, the behavior of DH is drastically different (Fig. 2).We then present reliable computations of length-scales associated with DH in an active system via three different methods (Fig. 3) and finally, illustrate the nontrivial aspects of activity (Fig. 4).In particular, we show that systems with the same relaxation time can have different DH; thus the dynamic length scale may not always control the relaxation dynamics in a glassy system, at least when energy is pumped at a microscopic length scale.
In this work, we study two different, but well-known, three-dimensional models of glassy systems: the Kob-Andersen binary mixture (3dKA) with 80 : 20 ratio of two types of particles and another binary mixture (3dHP) with 50 : 50 ratio of two types of particles with diameter ratio of 1.4.Results of 3dKA model are presented in this article and some of the results of 3dHP model are shown in SM, Fig. S10.Activity in our systems are introduced by designating 10% of the particles, chosen at random, as active.These active particles are then given constant forces of magnitude f 0 in random directions with the constraints that vector sum of these applied active forces adds up to zero; this ensures momentum conservation of the centre of mass in our simulations.The applied active forces are reshuffled every τ p time.In this work, we keep density and τ p = 1.0 constant and study the glassy properties as functions of temperature, T , and f 0 .We also extend active MCT, presented in Ref. [26], to study DH following the IMCT formalism presented in [47].More details of the models and the simulations can be found in the materials and methods section and in the SM.The details of the theoretical calculations will be presented elsewhere.

The relaxation dynamics
We first show that the relaxation dynamics of an active glass is similar to that of an equililibrium system at a suitably chosen T ef f .The overlap function, Q(t), and the mean-squared displacement (MSD), ∆r(t) 2 , are shown in Fig. 1(a) and in the inset respectively, as a function of time t (see Materials and Methods for definitions).At a particular T = 0.45, both Q(t) and MSD for the passive system (f 0 = 0) show the plateau, characteristic of glassy dynamics.As activity increases from 0, decay of Q(t) and growth of MSD become faster, as has also been found in earlier simulations [20][21][22] and experiments [16].The relaxation time, τ α , is defined as Q(τ α ) = 1/e.MCT for the active systems [26] predicts where γ is the same exponent as in a passive system, T M CT is the mode-coupling transition temperature in the absence of activity, and K, a constant.For the particular model of activity that we are using in this work, defined as SNTC in Ref. [26], we have K = Hτ p /(1 + Gτ p ) with H and G being two constants.We show the data of τ α (T ) for the equilibrium 3dKA model in Fig. 1 (b) (all simulation data, presented here are for this model, results for the 3dHP model are shown in SM Fig. S10).A fit of the MCT prediction, τ α ∼ |T − T M CT | −γ , with the simulation data gives γ ≃ 2.36 and T M CT ≃ 0.423.These numbers are similar to previously reported results for the passive system [48,49].Using these values, we fit one set of data for the active system and obtain K ≃ 0.015.We then plot Eq. 1 along with the simulation data in Fig. 1(c) and show that the theory agrees remarkably well with the simulation data.Note that once K is obtained from one set of data, there are no other free parameters in the theory.Moreover, Eq. 1 shows plot of Y = τ α (T − T M CT ) γ as a function of X = f 2 0 /(T − T M CT ) at different values of activity should follow a master curve: Y ∼ (1 + KX ) −γ .Fig. 1(d) shows excellent agreement between simulation data and the MCT prediction.
From Eq. (1), we also find that both the MCT transition temperature, T C , for the active system as well as the calorimetric glass transition temperature, T g , defined via τ α (T = T g ) = 10 6 , behave as T C,g = a − bf 2 0 , with a and b being two constants.As shown in SM Fig. S7, this prediction also agrees well with simulation data.The excellent agreement between the theory, Eq. 1, and simulation data with γ being the same exponent as for the passive system can also be interpreted as the relaxation dynamics of the active system being similar to that of the equilibrium glass at However, as we show below, an active glass behaves quite differently from an equilibrium glass.The nontrivial effects of activity manifest via the DH and the dynamic length scale, characterized via four-point correlation function, χ 4 (t).We now present the extension of active MCT to study DH via χ 4 (t) along with the simulation results for DH and dynamic length scale.

Active In-homogeneous Mode Coupling Theory (Active-IMCT)
To the best of our knowledge, MCT for an active system has not been extended to study the DH, the detailed derivation will be presented elsewhere, here we briefly outline the main result.We follow the IMCT formalism, developed in Ref. [47], and obtain the four-point correlation function in terms of a corresponding susceptibility via linear response theory.Since an active system is inherently out of equilibrium, one must study both the correlation and (integrated) response functions, Q(t) and F (t) respectively, as well as their corresponding susceptibilities.We impose a weak external field, ε, derive the equations of motion for Q(t) and F (t) in the presence of ε, and obtain the susceptibilities, χ Q (t) and χ F (t), defined as The equations of motions for χ Q (t) and χ F (t) are where, and where ∆(s) is the active force statistics: . κ is a system-dependent constant, set to κ = 2.0 in this work and a more detailed analysis will be presented elsewhere.As shown in Ref. [47], χ Q (t) ≡ χ 4 (t) is the desired four-point correlation function.These equations must be simoultaneously solved along with the equations for Q(t) and F (t), as presented in Ref. [26].

Dynamic heterogeneity in active glasses
We now present the simulation results showing the effect of activity on DH, characterized via χ 4 (t), defined in the Materials and Methods.Fig. 2(a) shows χ 4 (t) at T = 0.45 and different values of f 0 ; χ P 4 , the peak value of χ 4 (t), decreases with increasing f 0 .Similar behavior is also found within MCT as shown in Fig. 2(b) where we have used λ = 1.95 (the MCT transition is found at λ = 2.0 in our notation).Fig. 1(c) shows that f 0 decreases τ α and drives the system away from the glassy regime, therefore, decrease of χ P 4 with increasing f 0 is expected; but, the trend is different from equilibrium-like behavior.f 0 as well as the data for f 0 = 0.The lines are fit with a function χ P 4 ∼ (T − T C ) −ν .If the glassy dynamics of the active system resembles the equilibrium dynamics at a suitably defined T C , one expects ν to remain same at different activity and identical to that of the passive system.However, we find the exponent ν depends on activity: 0.92 for the passive system, whereas ν = 1.03 and 1.13 for f 0 = 2.0 and 2.5 respectively.Irrespective of the detailed values of ν, the important point of this analysis is that effect of activity on DH is quite different from the equilibrium-like dynamics.To probe this point from another perspective, we now study χ P 4 as a function of f 0 at different T (Fig. 2d).A fit of the function χ P 4 ∼ (1 + bf 2 0 ) −µ , where b is a constant, with the simulation data shows µ depends on T (inset, Fig. 2d).Thus, varying activity not-only changes T C , it also affects the exponents.The time τ peak , when χ 4 (t) attains its maximum, provides a measure of a relaxation time.As shown in the SM Fig. S8, τ peak is proportional to τ α .Thus, the relaxation dynamics, characterized via either of Q(t) or χ 4 (t), can be described in terms of a T ef f (Figs.1c and 1d).These results show quite distinctive effects of activity on relaxation dynamics and DH; whereas effect of activity can simply be treated in terms of a T ef f for the former, it cannot be done so for the latter.
Within MCT, the nontrivial effect of activity on the dynamic heterogeneity is contained in the last term in Eq. ( 7).This term is a direct consequence of the nonequilibrium nature of the system and absent in equilibrium.An important characteristic of activity is the giant number fluctuation (GNF) [36,37], where number fluctuation, δN , is related to the mean, N , as δN 2 ∼ N α with α > 1. GNF has been experimentally measured in Ref. [50] for the system of Ref. [45] and it has been argued that the nontrivial behavior of χ 4 (t) with activity might be related to GNF [46,50].The detailed mechanism of how GNF affects DH and why its effect on the relaxation dynamics remains equlibrium-like remains an intriguing question.The distinctive natures of the effects of activity on the relaxation dynamics and DH length scales imply they decouple from each other.Such a decoupling, even in equilibrium, has been suggested in [42], where it has been convincingly shown that the finite size effects of χ P 4 and τ α are not controlled by the same length scale.In fact, this is one of the first direct signatures of the existence of yet another length scale, the static length scale of amorphous order [51], related to structural relaxation time.It remains to be seen whether the strong decoupling of DH and structural relaxation in active glasses is an enhanced manifestation of the same equilibrium phenomena or distinct from it.

Dynamic Heterogeneity Length Scale
Having discussed the non-trivial effects of activity on the higher order correlation functions like χ 4 (t), we now turn our focus to compute the dynamic length scale, ξ D , associated with DH itself.Reliable estimate of ξ D in model active glassy systems is rare [52]; in this work, we compute and compare ξ D via three different methods -(i) finite size block scaling of χ 4 , (ii) block analysis of van Hove function and (iii) spatial correlation of displacements at structural relaxation time (see SM for further details): the agreement among them shows reliability of our results.We have used the method of block analysis that is developed to study the growing dynamical correlations in equilibrium supercooled liquids in [53,54].The simplicity and the novelty of the method makes it immediately useful for our current study as we need to compute ξ D not only as a function of T but also as a function of changing f 0 .Block analysis made it much easier for us to do these measurements for different values of f 0 with much less computational difficulty.Briefly, in block analysis method, a somewhat bigger system size with linear size L is simulated at the desired state points and then the whole system is divided into smaller blocks of length L B = L/n, where n is the number of blocks chosen, as schematically shown in Fig. 3(a).We typically choose largest blocks to have the linear size L/3, to ensure that other boundary effects do not influence the results.It is important to note that in this analysis, the blocks can be thought of in a Grand-Canonical ensemble, thus ensuring all possible fluctuations e.g.number of particles, density, and, T contribute in χ 4 (t) properly [53,55].Now, the dynamic susceptibility is computed for each such blocks as where N B is the number of blocks with size L B , n i is the number of particles in the i th block at time t = 0, and the window function w(x) = Θ(a − x) where Θ is the Heaviside step function and the value of the parameter a is chosen to remove the decorrelation arising from vibrations of particles inside the cages formed by their neighbours.The dynamic length scale is obtained from a detailed finite-size scaling (FSS) analysis of the dependence of the four-point dynamic susceptibility on the block size by assuming the following scaling ansatz, where, χ 0 (T ) = lim L B →∞ χ P 4 (L B , T ).The data for all temperatures can be collapsed to a master curve using the two parameters, χ P 4 (∞, T ) and ξ D .Fig. 3(b) shows χ P 4 (L B , T ) as a function of block length L B for various T for the passive system (f 0 = 0) and Fig. 3(c) shows the data collapse, Eq. ( 11), to obtain ξ D as shown in Fig. 3(h).In Fig. 3(d), we show similar data collapse for an active system with f 0 = 2.5.Similar analyses for all other values of f 0 are reported in SM along with other relevant discussion.The excellent quality of the data collapse observed in the analyses suggests the reliability of the extracted ξ D .However, in the absence of previous data to compare with, we reconfirm the reliability of the estimates of ξ D via two other methods as mentioned before.
It was demonstrated in Ref. [54] that ξ D can be obtained by analyzing the non-Gaussian behaviour of the self part of the van Hove function, G s (x, τ α ) (defined in Method section and elaborated in SM).In this method, one looks at the behaviour of the van Hove function as a function of spatial coarse-graining scale, L B , similar to the block analysis with the main difference, that here one computes the van Hove function of the displacement fields of the centre of mass (CoM) of the individual blocks.Intuitively, as the coarse-graining volume gradually increases, the system will eventually become spatially homogeneous once the scale of coarse-graining exceeds the typical dynamic heterogeneity length scale.Fig. 3(e) shows G s (x, τ α ) with increasing L B illustrating the change of nature of G s (x, τ α ), from non-Gaussian to Gaussian form, as L B increases.G s (x, τ α ) becomes Gaussian at a crossover scale L B which is different at each T and the T -dependence of this crossover scale has the same T -dependence of ξ D [54].To estimate this crossover length scale in an unbiased way, we computed the Binder cumulant of the distribution as shown in Fig. 3(f).It is clear that with decreasing T , the Binder cummulant increases strongly as L B decreases.As L B → ∞, Binder cumulant should reach 0 for at all T .We can do a one parameter scaling collapse (Fig. 3g) using the scaling ansatz B(L B , T ) = G (L B /ξ D ) where ξ D is the crossover length scale, same as the one used in collapsing the χ P 4 (L B , T ) data in the top panels.The quality of the data collapse shows that these two methods yield the same ξ D .
Finally, we have extracted the same length scale via yet another way used earlier in [56], by studying the spatial correlation of the displacement fields of the particles over structural relaxation time scale characterized by the correlation function g uu (r, τ α ), where r is the radial distance (see the definition in the Method section and further details in SM).In the inset of Fig. 3(h), the comparison of the length scales are shown.The excellent agreement again highlights the reliability of the obtained length scale.Fig. 3(h) shows the growth of ξ D for various activities as a function of temperature.Note the dramatic growth of heterogeneity length scale with increasing activity, which is one of the main results of this work and we are going to discuss the salient features of this observation in the next section.

Dramatic Effect of Activity on Dynamical Heterogeneity
We now discuss the dramatic effect of the non-equilibrium active force on the DH in an active glass.Within the DH picture of equilibrium glassy dynamics, it is beleived that ξ D controls the relaxation dynamics [1,38].We have already shown that activity affects DH in a way that is different from equilibrium-like behavior (Fig. 2c, 2d).Active-IMCT also suggests that, the memory kernel within MCT has a nontrivial effect of activity via the presence of the second term in Eq. (7).It is easier to demonstrate the nontrivial effects of activity by considering different systems whose parameters are chosen such that the typical relaxation times of the systems remain same.In an equilibrium system, DH, as it controls relaxation dynamics, should remain unchanged when τ α is same.However, the scenario for an active glass, surprisingly, turns out to be drastically different.
DH essentially refers to the coexistence of dynamic fast and slow moving regions in the same system; a set of adjuscent particles with similar τ α relaxes collectively and is known as cooperatively rearranging region (CRR) that The data shows that for a given τα, a system with higher f0 has larger ξD.(d) Fit of simulation data (symbols) for ξD with the power-law prediction of MCT (lines); the exponent γD is expected to be constant for an effective equilibrium-like scenario, however, we find that γD almost linearly increases with f0 (inset).(e) In simulation, we choose a set of T and f0 such that T ef f remains similar.(f) χ4(t) for the parameters as in (e); system with larger f0 has higher χ P 4 .Inset: All these systems have similar τα as seen from the plot of Q(t).(g) χ4(t), obtained from active-IMCT, also has similar behavior as in simulation, i.e., higher χ P 4 for larger f0 when Q(t), and hence τα, are similar (inset).(h) Variation of van Hove function for different activity; T is choosen such that the structural relaxation time is same.The increasing non-Gaussianity signifies the drastic increase in dynamic heterogeneity with increasing activity.(i) Increasing spatial correlation, as measured by the excess displacement-displacement correlation function, Γ(r, τα), with increasing activity.These correlation functions are again computed at the same relaxation time to highlight the growth of spatial heterogeneity with increasing activity (see text for details).
gives a measure of DH [57].Figs.4(a) and 4(b) show typical estimates of CRRs of the fast moving particles in an equilibrium system and an active system respectively, both having the same τ α = 10 3 (see SM, Sec. 5 for details).Whereas the CRR for the former is disconnected at this particular τ α , that for the latter becomes system-spanning.This qualitative comparison of typical CRRs in the two systems show greater DH in an active system at the same τ α .Now, for a more quantitative analysis, we show ξ D as a function of τ α in Fig. 4(c).As is evident, at a given τ α , a system with higher activity (that is larger f 0 ) has larger ξ D .Note the drastic increase of ξ D for the largest activity case, f 0 = 2.5, ξ D increases by a factor of ∼ 30 compared to its high temperature value, whereas in passive system, the growth is hardly by a factor of 4 − 5 in the comparable window of τ α .Our results are consistent with the large dynamic length scales found in simulations of a driven granular fluid [44].
MCT predicts power-law divergence of ξ D at T C , where T C is the mode-coupling transition temperature for an active system, obtained from fitting the data of τ α .Fig. 4(d) shows simulation data of ξ D as a function of (T − T C )/T C and the lines are fits with the function ξ D ∼ [(T − T C )/T C ] −γ D ; the data agree well with the power-law form in the range of T and f 0 where the simulations are carried out.An effective equilibrium-like description implies same γ D at different activity, as for the case of relaxation dynamics, Eq. (1).Consistent with the T -dependence of χ P 4 (Fig. 2d), we find γ D depends on activity and almost linearly increases with f 0 , as shown in the inset of Fig. 4(d); this result, again, highlights the nontrivial activity-dependence of ξ D .
To further illustrate the decoupling of DH and τ α in the presence of activity, we now show simulation data for Q(t) and χ 4 (t) at different T and f 0 in Fig. 4(f).The parameters are chosen in such a way that all the systems have very similar τ α , as is avident from the plots of Q(t) shown in the inset of Fig. 4(f).Analysis of τ α in terms of T ef f implies all the systems have similar T ef f , this is indeed true as shown in Fig. 4(e).If DH controls the relaxation dynamics, as in an equilibrium system, one expects χ 4 (t) for all the systems to overlap.However, as shown in Fig. 4(f), this is not the case; χ 4 (t) increases with increasing f 0 .Similar behavior is also found within our active-IMCT as shown in Fig. 4(g) with the inset showing Q(t) overlapping with each other for these chosen parameters.Thus, the effects of activity on the DH and relaxation dynamics are different and the unique relation between DH and τ α breaks down in an active glass.Similar conclusions can be drawn by studying the effect of activity on other dynamical quantities as well.Fig. 4(h) shows the increasing non-Gaussian behaviour in G s (x, τ α ) with increasing activity while the structural relaxation time is kept constant, indicating growing DH with increasing activity even when the structural relaxation time is same.Similarly, Fig. 4(i) once more demonstrates that spatial correlation, computed by the excess part of the displacement-displacement correlation function, Γ(r, τ α ), increases markedly with increasing activity, confirming the strong decoupling of DH and structural relaxation dynamics in active glasses.

CONCLUSIONS
To conclude, we have studied the dynamic heterogeneity in a glass-forming system in the presence of non-equilibrium active forcing.We find that the effect of activity on relaxation dynamics of an active glass can be described in terms of a T ef f , much like an equilibrium system.However, the notion that the effect of activity simply fluidizes the system with no other significant changes in its dynamical properties is not accurate and an active glass is fundamentally different from an equilibrium glassy system, the nontrivial effect of activity manifests via the DH.Through the intensive research of the last couple of decades, DH has emerged as one of the salient features of glassy systems [1,38,39] and is believed to control the relaxation dynamics.In the presence of activity, systems with the same τ α can have different DH and ξ D , thus, masking the unique relation between DH and the relaxation dynamics in an active glass and casting doubt on the presumed central role of DH in glassy dynamics.It is possible that the giant number fluctuations, which is a feature of active systems [17,36,37], affects DH and τ α differently.In any case, our results expose the incompleteness of our understanding on how DH affects relaxation dynamics and the fundamental properties of DH in glassy systems in general.
We have shown that four-point susceptibility that quantifies DH in a glassy system has a nontrivial behavior in an active glass.Considering glass transition as a critical phenomenon, activity-dependence of the exponents of the power-law behavior of both χ P 4 and ξ D posits a serious challenge in developing a microscopic theory for such a system.Yet, the active-IMCT does predict the nontrivial effects of activity on DH, which, the theory suggests, comes from an activity-specific term within the memory kernel, Eq. ( 7).To the best of our knowledge, exploration of the effect of GNF on the glassy dynamics is scant [50].Studying GNF in a dense glassy system may not be straightforward, however, the vestige of the effect should be there in the system.Due to its relevance in several biological systems, active glass is an exciting system to investigate this effect as an interesting future direction.
An important result of the current work is the dramatic growth of DH in an active system, consistent with simulation results [44] on driven granular fluids, which is a prototype active system [11,17,36].In particular, we show that dynamic heterogeneity grows dramatically with increasing activity even when τ α remains constant.Recent experiments on confluent cellular monolayer have found increased χ P 4 via the application of GTPase RAB5A, even when the system fluidizes, i.e., τ α decreases [45,46].Application of RAB5A may affect a confluent monolayer in several ways, such as affecting the junctional proteins, and thus, changing inter-cellular interaction [58] as well as affecting the motor proteins [59].Via a combination of theory and experiment, it was convincingly shown in Ref. [50] that at least one important effect of RAB5A on the monolayer is increasing cellular motility.Thus, these exciting results, though puzzling from equilibrium glassy perspective where increased χ P 4 implies higher τ α , are entirely consistent in the light of our theory.
Finally, we want to highlight that this work for the first time demonstrates that the dynamical heterogeneity in active glassy systems is inherently different from their equilibrium behaviour.Via a combination of theory and simulation, we have shown that whereas the relaxation dynamics can be described as an equilibrium-like behavior at a certain T ef f , the DH has completely different behavior.This also highlights the decoupling of relaxation dynamics and DH, an effect that has also been stressed in some equilibrium systems [42].Studying the active glass, thus, in general promises deeper understanding of the role of DH in glassy dynamics.

MATERIALS AND METHODS
We study two different glass-forming models, one is Kob-Andersen binary mixture with 80% A-type and 20% B-type particles interacting via the Lennard-Jones pair potential where r is the distance between two particles and the indices i and j can be A or B. The values of σ ij and ǫ ij are chosen to be: σ AB = 0.8σ AA , σ BB = 0.88σ AA , ǫ AB = 1.5ǫAA, and ǫ AB = 0.5ǫ AA .We set a cutoff in the potential at r ij = 2.5σ ij and shift it accordingly.We set the unit of length and energy as σ AA = 1 and ǫ AA = 1 and fix the overall density ρ at 1.2.We use a quadratic polynomial to make the potential and its first two derivatives smooth at the cutoff distance.We chose another model that interpolates between finite-temperature glasses and hard-sphere glasses and has been studied extensively in the context of jamming physics.This is a 50 : 50 binary mixture with diameter ratio of 1.4.Details of the potential for this model is given in the SM.
We have done all our simulations in constant number of particle (N ), volume (V ) and temperature (T ) (NVT) ensemble.To keep the temperature constant we use Gaussian operator-splitting thermostat [60] throughout our simulation.Following Ref. [21], we introduce activity by making 10% of the particles active (see Sm for details).In this work, we keep τ p = 1.0 fixed and study the effect of activity as a function of f 0 alone.The opposite trends of activity, discussed in Ref. [61], is not important in our study since both these models have same trends as a function of f 0 .

DYNAMICAL QUANTITIES
To characterise the dynamics of glassy systems under various active forcing, f 0 , we measure the following dynamical quantities.Mean-squared displacement: The mean-squared displacement (MSD) ∆r(t) 2 is defined as Overlap correlation function: The two point overlap correlation function Q(t) is defined as where ... represents the average over time origin t 0 as well as 10 independently generated configurations, N is the total number of particles in the system and the parameter a is associated with the typical vibrational amplitude of the caged particles.Throughout our analysis, we have used the parameter a = 0.3 and we have verified that our results do not change with moderate variation in a. Four-point dynamical susceptibility: The four-point dynamical susceptibility χ 4 (t) is defined in terms of the fluctuations of the two point overlap correlation function as Excess displacement correlation function: The excess displacement correlation Γ(r, ∆t) defined as Γ(r, ∆t) = g uu (r, ∆t) where g(r) is the radial pair correlation function, defined as g(r) = 1 ρN N i,j=1,j =i δ(r + r i (0) − r j (0)) , and the spatial correlation of the particle displacements, g uu (r, ∆t), defined as g uu (r, ∆t) = N i,j=1,j =i u i (t, ∆t)u j (t, ∆t)δ(r− | r ij (t) |) 4πr 2 ∆rN ρ u(∆t) where u i (t, ∆t) =| r i (t + ∆t) − r i (t) | is the scalar displacement of the particle between time t and t + ∆t.
van Hove correlation function: The self part of the van Hove correlation function is defined as Author Contributions: SK designed the research project.KP performed the research.SKN developed the active MCT formalism.KP, SKN and SK analysed the data.SK and SKN wrote the paper.

FIG. 1 .
FIG. 1. Characterization of relaxation dynamics.(a) Decay of the overlap function, Q(t), and growth of MSD, ∆r(t) 2 , (inset) becomes faster with incresing activity.(b) Symbols are τα as a function of T for a passive system and the line is a fit with the MCT prediction τα ∼ (T − TMCT ) −γ that gives TMCT = 0.423 and γ = 2.36.(c) τα as a function of f0 at different T , symbols are simulation data and the lines are active MCT prediction (not individual fits).(d) Equation (1) predicts a master curve when τα(T − TMCT ) γ is plotted as a function of f 2 0 /(T − TMCT ) at different T and f0: the simulation data (symbols) agrees well with the theoretical prediction (line).

FIG. 2 .
FIG. 2. Characterization of χ4(t) as a function of T and activity.(a) Simulation data for χ4(t) at different f0 and T = 0.45.χ P 4 , the peak value of χ4(t), decreases with increasing activity.(b) Active-IMCT results for χ4(t) at different f0 and λ = 1.95; theory shows similar results as in simulation.(c) Symbols are simulation data for χ P 4 as a function of T at different f0 and lines are fits with a power-law where the exponent depends on f0 (see text).(d) χ P 4 as a function of f0 at different T ; fitting the simulation data (symbols) with a power-law form gives the exponent µ that depends on T (inset).

FIG. 3 .
FIG.3.Estimation of dynamic heterogeneity length scale.(a) Schematic of block analysis method showing how the whole system is subdivided into smaller blocks to do the finite size scaling analysis.(b) The block size dependence of peak height of χ4(t), χ P 4 , for the passive system.(c) The finite size scaling collapse of the data in panel (b) to obtain the dynamic heterogeneity length scale, ξD (Eq.11).(d) Similar scaling collapse, as in (c), of χ P 4 for an active system with f0 = 2.5.(e) Block size dependence of the self part of the van Hove correlation function, Gs(x, τα), computed at structural relaxation time.At small LB, Gs(x, τα) shows non-Gaussian behaviour and becomes progressively Gaussian with increasing LB.Gs(x, τα) becomes Gaussian when LB exceeds a crossover length scale that gives ξD.(f) Binder cummulant of the block averaged Gs(x, τα) at different T .(g) Scaling collapse of the Binder cummulant computed in panel (f) to obtain ξD.(h) ξD as a function of scaled temperature, T /TK , where TK is the Kauzmann temperature at various activities.Note that the same length scale collapses both the χ P 4 and Binder cummulant of van Hove function data.Inset: Comparison of different estimates of ξD, obtained via block analysis of χ P 4 and using spatial displacements correlation function, guu(r, τα) (see text for details).The agreement is very good suggesting the robustness of the methods as well as the estimate of ξD.

FIG. 4 .
FIG.4.Nontrivial effects of activity on glassy dynamics.(a-b) Qualitative measures of CRR showing the fraction of fast relaxing particles, CRR provides a typical estimate of DH, for a passive (a) and an active system (b) at same τα = 10 3 and shows larger DH in an active system.(c) ξD as a function of τα at different f0.The data shows that for a given τα, a system with higher f0 has larger ξD.(d) Fit of simulation data (symbols) for ξD with the power-law prediction of MCT (lines); the exponent γD is expected to be constant for an effective equilibrium-like scenario, however, we find that γD almost linearly increases with f0 (inset).(e) In simulation, we choose a set of T and f0 such that T ef f remains similar.(f) χ4(t) for the parameters as in (e); system with larger f0 has higher χ P 4 .Inset: All these systems have similar τα as seen from the plot of Q(t).(g) χ4(t), obtained from active-IMCT, also has similar behavior as in simulation, i.e., higher χ P 4 for larger f0 when Q(t), and hence τα, are similar (inset).(h) Variation of van Hove function for different activity; T is choosen such that the structural relaxation time is same.The increasing non-Gaussianity signifies the drastic increase in dynamic heterogeneity with increasing activity.(i) Increasing spatial correlation, as measured by the excess displacement-displacement correlation function, Γ(r, τα), with increasing activity.These correlation functions are again computed at the same relaxation time to highlight the growth of spatial heterogeneity with increasing activity (see text for details).