From motility-induced phase-separation to glassiness in dense active matter

Dense active systems are widespread in nature, examples range from bacterial colonies to biological tissues. Dense clusters of active particles can be obtained by increasing the packing fraction of the system or taking advantage of a peculiar phenomenon named motility-induced phase separation (MIPS). In this work, we explore the phase diagram of a two-dimensional model of active glass and show that disordered active materials develop a rich collective behaviour encompassing both MIPS and glassiness. We find that, although the glassy state is almost indistinguishable from that of equilibrium glasses, the mechanisms leading to its fluidization do not have any equilibrium counterpart. Our results can be rationalized in terms of a crossover between a low-activity regime, where glassy dynamics is controlled by an effective temperature, and a high-activity regime, which drives the system towards MIPS. Emergent behavior of dense living materials shows coexistence of both equilibrium and non-equilibrium features, such as glassy states and motility-induced phase separation. Here, the authors study transition between these two phases in a model of dense, disordered, soft active materials and find that the mechanism leading to fluidization from the glassy phase do not have an equilibrium counterpart.

A ctive self-propelled particles have been often employed as a starting and minimal model for capturing collective behaviors in biological and living materials, examples run from biological tissues to dense drops of ants [1][2][3][4][5] . Moreover, active agents might show a high degree of heterogeneities in their microscopic characteristics inside a given population. However, in modeling biological tissues, usually, it is assumed that each cell has the same mechanical and geometrical properties 2 . Cellular differences play a crucial role in many biological processes as in the case of cancer development, where cell heterogeneity is a common feature 6 . Phenotypic heterogeneity has functional consequences that impact the ability of microbes to adapt to different environments 7 and mechanical heterogeneity changes the collective properties in simple models of biological tissues 8 .
Living systems composed of a collection of self-propelled agents are thus naturally characterized by a certain degree of variability in size, shape, mechanical properties, and motility parameters. Often, in the system of interest, only average values are available: it is thus reasonable to introduce suitable probability distributions for those quantities.
In this work, we will focus our attention on the effect of quenched fluctuations in the typical active agent size. To take into account this aspect, we perform numerical simulations of Active Brownian disks of heterogeneous size. To make a contact with what is known on equilibrium systems, we study the phase diagram of repulsive Active Brownian Particles (ABP) composed of a continuous, i.e., polydisperse, mixture 9 . Systems of polydisperse self-propelled particles play an important role for our understanding of dense and active materials [10][11][12][13][14][15][16] . On the other hand, understanding the glass transition in active matter requires a huge theoretical and numerical effort and many progresses in the framework of Mode-Coupling Theories have been done during the last years [17][18][19][20][21] .
The model system introduced here is particularly suitable for studying the interplay between Motility-Induced Phase Separation (MIPS) and glassy dynamics [22][23][24][25][26][27] . MIPS is a well-established phenomenon in Active Matter whose main features are very robust against different microscopic active dynamics 28,29 and dimensionality of the space 30 . We shall see that glassiness occurs only at low activity, although the system in the MIPS phase arranges in dense patches that do not provide any evidence of positional order (and thus glassy behavior might develop).
Using as a control parameter the persistence length ℓ and the packing fraction ϕ, we study the phase diagram of polydisperse ABP and investigate the properties of the system at high packing fractions, where we document a glass transition driven by the persistence length. The latter is signaled by a dynamical slowing down of the time correlation function of the hexatic order parameter, with a decoupling between density and hexatic fluctuations. We show that, in the glassy region, there is a coexistence of hexatic and liquid phases. We then look at the statistical properties of the inherent structures obtained from stationary configurations at different parent activities. In this way, we can provide an estimate of the distance, computed in a probabilistic sense, between a typical steady-state configuration, that results from the non-equilibrium dynamics in the presence of self-propulsion, and the corresponding configuration that minimizes the mechanical energy. We show that glassy configurations are almost inherent configurations. Comparing dense configurations obtained for higher persistence length, we document a discontinuous crossover approaching MIPS.
Finally, we explore the concept of effective temperature in connection with the vibrational density of states of instantaneous configurations. In general, the eigenvalues of the dynamical matrix can be either positive or negative in the case of instantaneous configurations. Since negative eigenvalues are connected to unstable directions in the instantaneous energy landscape, they can provide information on barrier crossing [31][32][33] .
We show that the magnitude of the largest negative eigenvalues of the dynamical matrix, i.e., the "faster" unstable mode of the system, encodes information about the break down of the concept of effective temperature in active glasses. Those unstable directions can be surfed by the active system for removing geometrical frustration and bringing the system into a fluid phase. This result suggests that, even though glassy configurations at low activity are basically indistinguishable by their equilibrium counterpart, the mechanism of fluidization in active systems might be distinct from that in equilibrium glasses, where mean-field models point the attention out on marginal directions in the potential energy landscape rather than the unstable ones 34 .
In the present work, as main results we obtain that (i) the MIPS region extends up to very large packing fractions, (ii) there is no hint of dynamical slowing down in the MIPS regime, even if the dense MIPS phase is an assembly of polydisperse particles at high packing fraction, (iii) the mechanisms of fluidization of the active systems might be quantitatively connected with unstable modes of the instantaneous spectrum.

Results
Absence of structural arrest at high activity. We start our study considering a system composed of fully polydisperse Active Brownian Particles (ABPs, details about the microscopic model are provided in Materials and Methods). A qualitative view of the phase diagram is provided in Fig. 1 where we report the snapshots of typical stationary configurations. As one can see, the system undergoes MIPS and separates into a dense and dilute phase for large enough values of the persistence length ℓ. We observe that the MIPS regime extends until the largest packing fraction simulated, i.e., ϕ = 0.78. As we shall study in the next sections, the high density regime is characterized by a glassy phase at small persistence lengths.
We start our quantitative analysis by looking at the statistical properties of the relevant coarse-grained fields, i.e., the local packing fraction field ϕ(x, y, t), which measures the packing fraction around the point of coordinate (x, y) at time t, and hexatic order parameter ψ 6 (x, y, t), which provides a local measure of the orientational order in (x, y) at time t. Details about their definition and their computation are provided in the section Methods. For monitoring the presence of MIPS, we look at the behavior of the probability distribution function of ϕ(x, y, t), PðϕÞ ¼ hδðϕ À ϕðx; y; tÞÞi t (see Materials and Methods for details). To investigate the emergence of a disordered glassy phase, we study the relaxation time of density fluctuations through the intermediate scattering function F self (q, t) and the time-correlation of the hexatic order parameter C ψ (t) (their definition is provided in Methods). The resulting phase diagram is shown in Fig. 2a where the color map indicates the structural relaxation time measured through C ψ (using τ α defined from the decay of F self does not introduce any qualitative changing in the phase diagram). As one can appreciate, although in the MIPS region the system separates into a dilute and a dense phase (with a dense phase reaching packing very high packing fractions, as shown by the behavior of PðϕÞ), the glassy region is limited to high packing fraction and small persistence length. It is worth noting that, in systems composed of Active-Ornstein Uhlenbeck particles, where the natural control parameters are the correlation time of the noise and its strength 35 , keeping constant the diffusion constant, one can obtain a glassy regime that goes to very large persistence lengths 16 .
The behavior of PðϕÞ crossing the MIPS region is shown Fig. 2b. The distribution develops the typical double-peaked structure due to the presence of two coexisting phases. Since the potential is very steep and soft, the tail of the peak at high packing fractions can reach ϕ ≳ 1. The same is true even at very high densities, as documented in Fig. 2c. The qualitative features of the phase diagram are consistent with those observed in athermal self-propelled disks 17 . We checked that the system remains always in disordered, liquid-like configurations, in the whole range of parameters, by looking at the radial distribution function g(r), as shown in Fig. 2d. g(r) does not reveal any hint of crystallization, even at the largest packing fraction ϕ = 0.78.
Polydisperse disks introduce geometrical frustration in the system, inhibiting crystallization. Although at high density and small persistence length a glassy regime develops, heterogeneous regions where the system develops hexatic patterns survive over short length scales (of the order of a few particle radii) at high densities. The presence of these heterogeneous zones of high and low ψ 6 values determine a region of the phase diagram that is delimited by the dashed blue line in Fig. 2a, and that we identify with the label Liquid/Hexatic. To investigating this feature, we explore the presence of hexatic patches, typical of twodimensional systems 28,[36][37][38] , by looking at the distribution Pðjψ 6 jÞ ¼ hδðjψ 6 j À jψ 6 ðx; y; tÞjÞi t , shown in Fig. 2e (ℓ = 0.1). At high densities, the distribution becomes double peaked indicating the presence of hexatic patterns.
Crossing the dashed blue line in Fig. 2a, Pðjψ 6 Þj changes from single to double-peaked. The Fluid/Hexatic transition seems to develop a reentrance in the phase diagram, as suggested by the behavior of the orientational order parameter 〈|ψ 6 |〉 all along the phase diagram (see Supplementary Note 1, where the contour plot of 〈|ψ 6 |〉 is shown). This is highlighted in Fig. 2f where we report the height of the second peak, i.e., P M max jψ 6 j > 0:5 Pðjψ 6 jÞ, as a function of ℓ at ϕ = 0.79. Within the statistics considered here, P M is a non-monotonous function of ℓ, i.e., P M decreases as ℓ increases towards the fluid phase and, eventually, develops a bump in the MIPS region. This behavior suggests that hexatic patches might be more pronounced in the glassy and MIPS region than in the liquid region.
Stability of MIPS against quenched fluctuations. In the previous section we have introduced the phase diagram of the fully polydisperse model. Now we are going to investigate the impact of geometrical frustration on the MIPS dense phase. In order to quantify the effect of particle heterogeneities on MIPS, we have performed numerical simulations deep in the phase-separated region for a larger system size, i.e., ℓ = 100, ϕ = 0.44, L = 120〈σ〉, and by varying the fraction f ∈ [0, 1] of polydisperse disks. This means that at f = 0, the system is monodisperse, i.e., σ i = 〈σ〉, ∀ i.
Typical stationary configurations are shown in Fig. 3. A quantitative analysis shows that, although Pðjψ 6 jÞ remains double-peaked at any f value, a small fraction of heterogeneous disks causes a huge decrease of the height of the second peak (see Fig. 4a). This fact can be quantified further by looking at P M defined before. As one can see in the inset of Fig. 4a, P M dramatically decreases as f increases. The reduction of orientational order can be appreciated by the study of g 6 (r). The typical behavior of g 6 (r) as f increases from 0 to 1 is presented in Fig. 4b. Although, because of finite size effects, it is hard to identify a clear power-law behavior at f = 0 39 . As f increases, g 6 (r) tends to decay exponentially, i.e., g 6 (r)~e −r/ξ , with a correlation length of the order of a few particle sizes that tends to screen any power-law decay. The fact that the system loses orientational order on large scales does not ensure the absence of a dynamical arrested phase. The latter is excluded by the behavior of the relaxation time of C ψ (the color map in the phase diagram for f = 1 shown in Fig. 2a) that provides clear evidence of a liquid-like dense phase (excluding any kind of glassy dynamics in the dense phase of MIPS). However, such structural change in the hexatic patches is not accompanied by any dramatic change on the density demixing due to MIPS. This is shown in Fig. 4c, where we report PðϕÞ for different values of f. We observe that height of the highdensity peak decreases as the fraction of heterogeneous disks increases. However, the position of the two peaks remains much more stable, indicating that both, the dense and the dilute phase, remain at almost the same average densities, i.e. the MIPS coexistence region remains largely unaffected by polydispersity. This effect can be quantified looking at height of the high-density peak of PðϕÞ, as it is shown in the inset of Fig. 4c. Although the peak slightly reduces its height as f increases, the distribution develops a fat tail for high ϕ values indicating a robust phase separation in dense and dilute regions. Particle colors indicate their radius (increasing values from violet to yellow). For large persistence lengths, the system undergoes Motility-Induced phase separation (MIPS) (snapshots with green background). The MIPS region extends to large packing fractions, i.e., ϕ ≃ 0.8. In the dense regime, the system behaves as a glass at small persistence lengths, then melts to a fluid state, and finally undergoes MIPS.
For gaining insight into the nature of the dense phase, we have computed the inherent structures at different parent activities. In practice, we consider an instantaneous configuration taken in the stationary state at ϕ = 0.79 for a given value of ℓ. We then compute the corresponding inherent configuration by minimizing the configurational energy Φ that is solely due to the pair interactions.
We start our study comparing the structural properties of inherent and instantaneous configurations. In order to perform a quantitative analysis, we compared the probability distribution functions PðϕÞ and Pðjψ 6 jÞ resulting from the inherent structures, i.e., P IS , with those obtained from instantaneous configurations, i.e., P inst .
The behavior of P IS;inst ðϕÞ is shown in Fig. 4d for ℓ = 0.001. In Fig. 4e the same observables are evaluated for ℓ = 200. As one can see, inherent and instantaneous configurations are quite similar at low persistence length (ℓ = 0.001), proving that the system reached a mechanically stable configuration that is slightly perturbed by activity. Different is the situation at high parent activities, i.e., ℓ = 200, Fig. 4e, where instantaneous configurations taken in the stationary state reveal the presence of MIPS through a broad tail in P inst ðϕÞ. On the other hand, the inherent configurations are homoheneous and thus they produce a Gaussian distribution centered around ϕ = 0.79.
For making quantitative progresses we measure the distance between P IS and P ins using the Kullback-Leibler divergence D KL ½P IS jP inst (see Methods). The result is shown in Fig. 4f for both φ and ψ 6 . We obtain that D KL ½P IS ðφÞjP inst ðφÞ is small in the glassy state, at small ℓ, indicating that instantaneous and inherent configurations are almost identical, and even decreases as ℓ  increases. However, as the system approaches MIPS, D KL jumps to higher values. This is because inherent configurations are homogeneous while instantaneous ones tend to decompose into two phases because due to MIPS. Less informative is D KL ½P IS ðjψ 6 jÞjP inst ðjψ 6 jÞ that maintains small values showing a smooth crossover on intermediate persistence length, i.e., ℓ~0.1. This crossover reveals measurable differences in the orientational order between inherent and instantaneous configurations.
Activity-driven dynamical arrest. We now explore the high packing fraction region ϕ = 0.79 in the case of fully polydispersisity, i.e., f = 1. We study the behavior of the self-part of the intermediate scattering function F self (q, t), and the time-correlation function of the hexatic order parameter C ψ (t). We monitor F self (q max , t), where q max is the wave vector of the first peak of the static structure factor, for increasing values of ℓ, see Fig. 5a. As one can see, activity tends to melt the system. However, as well as in the case of other two dimensional glassy systems 37 , F self (q, t) does not show a clear two-step decay, indicating that the cage effect is not as strong as in the case of a three-dimensional system. Looking at the sample-to-sample fluctuations of F self (q, t), we obtain a dynamical susceptibility χ 4 (q, t) that behaves as in the case of an ordinary supercooled liquid (Fig. 5f) 40 . The broad peak in χ 4 (q, t) reveals the presence of dynamical heterogeneity, as shown Fig. 5d, e, for ℓ = 0.001, 0.01, respectively, where we show the map of displacements, i.e., each arrow indicates the typical displacement performed by a particle on the time scale τ χ , with τ χ obtained from the peak of χ 4 (q, t). It is worth noting that τ χ behaves quantitatively in the same way as τ α , i.e., the structural relaxation time obtained from F self , (see Fig. 5c).
A more clear two-step decay is visible in C ψ (t) (Fig. 5b). This fact suggests the slowest structural degree of freedom of the system is the hexatic order parameter rather than the density. This claim is supported by the evidence that there is a bifurcation of the relaxation time scales of the two observables approaching the glass transition (Fig. 5c). We notice that similar differences in relaxation of F self and C ψ are documented in two-dimensional equilibrium glasses 41,42 where Mermin-Wagner excitations play an important role 43 . Because of the decoupling between τ χ and τ ψ , the order parameter ψ 6 remains frozen on the typical time scale of dynamical heterogeneities, as highlighted by the color map in (d) and (e).
It is worth noting that, as it has been shown in Fig. 2f, the height of the peak at higher |ψ 6 | values tend to reduce is height as the persistence length ℓ increases reaching a minimum at ℓ~1. As we will discuss later, this characteristic value of ℓ is related with the properties of the instantaneous normal modes. The same scenario has been observed at lower densities, i.e., ϕ = 0.66.
Instantaneous normal modes and the breaking of the effective temperature. The concept of effective temperature helps to rationalize some features of active systems [44][45][46][47][48][49][50][51][52] . In colloidal glasses driven by thermal noise, because of the caging effect, particles spend most of the time vibrating around their equilibrium position until cooperative rearrangements allow the particle to escape from the cage 53 . This equilibrium-like picture can be reasonably employed also in the case of an active glass whenever the active motion causes only vibrations and thus rattling inside the cage, i.e., in the small persistence length regime ℓ ≤ 〈σ〉. In order to make progress, we consider harmonic vibrations around an equilibrium inherent configuration. The harmonic vibrations define a set of eigenfrequencies ω 2 ν ¼ λ ν , with λ ν the the ν − th eigenvalue of the dynamical matrix 46,49,50 . We can thus introduce the following (mode-dependent) effective temperature 49 (see Supplementary Note 2) 35,54,55 ) T eff ðλ ν ; 'Þ ¼ v 2 0 2ðv 0 ='þμλ ν Þ , with μ being the particles' mobility. At variance with early studies, here we are going to provide an alternative interpretation of T eff (λ ν , ℓ) that allows to gain insight into the active glass transition.
Let r Ã ðr Ã 1 ; r Ã 2 ; :::; r Ã N Þ be the inherent configuration of the system so that ∂ i Φ r Ã ¼ 0. The stability of such a configuration is written in the eigenvalues λ ν of the Hessian matrix H ij ¼ ∂ i ∂ j Φ. If all the eigenvalues λ ν are strictly positive, the configuration lays in a minimum of the potential energy landscape. Let us introduce the persistent time τ = ℓ/v 0 . As soon as the system develops instabilities, i.e., directions with negative curvature in the energy landscape, there will be a critical value τ c that causes a divergence in T eff (for a given negative value λ ν , we can write λ ν = − |λ ν | and T eff loses completely sense for τ − μ|λ ν | = 0) so that, looking at instantaneous configurations that naturally develop unstable directions, T eff becomes ill-defined for increasing values of τ. This critical value corresponds to the largest negative eigenvalue λ ðÀÞ max , i.e., τ c ¼ jλ ðÀÞ max j À1 . For typical configurations in the liquid regime, the energy spectrum always develops negative eigenvalues, making the concept of effective temperature meaningful only in the limit τ ≪ τ c .
Replacing ℓ in favor of τ, we thus stress that not only T eff (λ ν , τ) makes sense for τ → 0, but also that the breaking of the concept of effective temperature is bounded to the structure of the energy landscape. For checking this fact, we have computed the spectrum of instantaneous and inherent normal modes. The corresponding density of states D inst;IS ðωÞ are shown in Fig. 6a. As a standard procedure, imaginary eigenmodes λ (−) have been mapped to negative frequencies, i.e, ω ðÀÞ ¼ À ffiffiffiffiffiffiffiffiffiffi ffi jλ ðÀÞ j q . At variance with instantaneous configurations, the frequency spectrum D IS ðωÞ of inherent configurations does not contain any negative frequency, meaning that it is linearly stable. As ℓ increases, the average number of negative frequencies 〈n (−) 〉 increases (Fig. 6b). Both 〈n (−) 〉 and the (slowest) structural relaxation time τ ψ approach a plateau value at the same crossover value ℓ c~1 .
For connecting the crossover with instabilities in the energy landscape, we look at the behavior of the largest negative frequency ω ðÀÞ max as a function of ℓ, as shown in Fig. 6c. The estimate from the effective temperature breaking is τ c jλ ðÀÞ max j ' 1 as one can see, we correctly obtain ℓ c~1 . We repeated the analysis at ϕ = 0.66 and the result fairly match the prediction (see Figs. 6c and 2b). As in the case of equilibrium glasses, we can thus relate the features of the energy landscape with dynamical slowing down 34,[56][57][58][59] . However, at variance with equilibrium, our picture suggests that the fluidization of glassy configuration in active system is connected to negative directions in the energy landscape rather than the marginal ones.

Discussion
Living systems as bacterial colonies or biological tissues are typically composed of self-propelled units of different sizes and different self-propulsive forces. In this work, we studied how this quenched disorder impacts the collective behavior of active particles. As a model system, we have considered a mixture of Active Brownian disks of different radii. We have explored the phase diagram using as control parameters the persistence length ℓ of the active motion and the packing fraction ϕ. Similar to the case of monodisperse active systems, we have documented the presence of a MIPS coexistence region for large persistence lengths. The phase separation extends from relatively small packing fractions (ϕ~0.2) to very high densities ϕ~0.8. Geometrical frustration due to geometrical heterogeneities makes the system considered here a good model of glass 9 . We showed that this peculiarity did not get lost because of activity. In particular, we documented a glass transition at high density characterized by decoupling between the relaxation time of density fluctuations and that of the hexatic order parameter.
The study of the inherent structures has unveiled a deep connection between stationary configurations at small persistence lengths and the corresponding inherent ones. This connection turns out to be lost as ℓ increases. This is because MIPS, although it might be interpreted in terms of effective potential and equilibrium-like spinodal decomposition, it is an intrinsic nonequilibrium phenomenon caused by the self-propelled motion which cannot be completely reduced to an effective equilibrium description.
The model introduced here is particularly suitable for studying the consequence of the competition between glassy dynamics, peculiar of equilibrium dynamics and known to play a pivotal role in dense active systems 60,61 , and MIPS, a typical feature of active systems whose importance in bacterial colonies has been recently documented 62 .
In this work, we have developed a deeper understanding of the concept of effective temperature in active matter, based on the analysis of the spectrum of vibrations around its inherent structures. Using a simple picture based on the definition of effective temperature through an extension of the energy equipartition theorem 45 that can be applied to dense active systems 46,49,50 , we found that the breaking of the concept of effective temperature in active matter can be linked to the properties of the energy landscape of the instantaneous configurations. Since an instantaneous configuration is not optimal, the corresponding vibrational spectrum contains always negative eigenvalues. We observed that the magnitude of the largest of them can be used as a proxy, i.e., τ c $ 1=λ ðÀÞ max , for estimating the crossover between active liquid, from the point of view of the relaxation time, and a active glassy regime, i.e., where the structural relaxation time starts increasing dramatically. The estimated value τ c~1 is also in agreement with a heuristic argument that estimates the validity of the concept of T eff for vibrations smaller than the typical particles size τ ≪ 〈σ〉/v 0~1 . Our approach, which is based on the properties of the vibrational energy landscape, provides an interpretation of the mechanisms leading to the fluidization of frustrated glassy configurations. In particular, our picture suggests that, in spite the similarities between glassy dynamics in active and equilibrium-driven systems, the fluidization of the active glassy state is due to microscopic mechanisms different from those in equilibrium, since the important parameter might be the magnitude of the largest negative eigenvalue rather than the number of quasi-stationary points 63 . As a future direction, it might be interesting to study other aspects of the topological crossover in active glasses 64 and how it is connected with nontrivial spatial correlations of the velocity field 65,66 . Moreover, while we have observed that the positional order in MIPS is strongly modified by the presence of polydisperse disks, less clear is how the phase boundaries of the MIPS region change with varying f.
Finally, the lack of any dynamical slowing down in the MIPS phase is consistent with the breaking of the concept of effective temperature. In particular, the MIPS regime lives always at high persistence lengths, and thus high activity, i.e., in a regime where non-equilibrium effects are important. This fact has strong implications, for instance, it follows that density and geometrical frustrations are not the only minimal ingredients for driving a dynamical arrest in Active Matter. More in general, MIPS seems to be incompatible with dynamically arrested phases. It seems Fig. 6 Instantaneous normal modes of fully polydisperse Active Brownian particles. a DðωÞ of inherent configurations (blue curve) and instantaneous configurations at different persistence length (ℓ = 10 −4 , 1, red and green, respectively.) b The average number of negative frequencies (red circles, error bar smaller than symbols, average over 100 independent configurations) approaches a plateau value at large ℓ and undergoes a crossover for ℓ~1 at the same point where the structural relaxation time (green squares) starts to grow for decreasing values of ℓ (blue-shaded area). c The largest negative frequency in unit of persistence time τ (red circles) as a function of persistence length ℓ is order 1 at the crossover. Magenta symbols refer to lower packing fraction, i.e., ϕ = 0.66. Green squares are the structural relaxation time. The dashed blue line is τjω À max j 2 ¼ ' (error bar indicates the statistical uncertainty over 100 independent configurations). crucial to understand whether or not the crossover from equilibrium-like glass and active regime might be captured by a Mode-Coupling Theory approach that suggests small vanishing values of the self-propulsion velocity of active tracers near the glass transition 18 .

Methods
We perform numerical simulations of a two dimensional Active Brownian particles composed of N disks of different sizes, interacting via a short-range power-law potential, and confined on a square box of length L with periodic boundary conditions.
Microscopic model. Due to their robustness against crystallization, polydisperse mixtures result to be good candidates as glass formers in a wide range of temperatures, even below the dynamical arrest temperature 9,67-71 . Here we consider a polydisperse mixture in two spatial dimensions where Active Brownian disks of different diameters σ i interact through a pair potential with r ij ≡ |r i − r j |. We set the softness exponent to n = 12. The coefficients c 0 , c 2 , and c 4 are chosen in a way that vðr c Þ ¼ v 0 ðr c Þ ¼ v 00 ðr c Þ ¼ 0, where we have introduced the notation v 0 ðrÞ ¼ dv dr and v 00 ðrÞ ¼ d 2 v dr 2 . For suppressing the tendency to demix, we consider non-additive diameters σ ij ¼ 1 where ϵ tunes the degree of non-additivity 9,72 . The cutoff is r c = 1.25σ ij , and ϵ = 0.2. The particle diameters σ i are drawn from a power law distribution P(σ) with hσi R σ max σ min dσ PðσÞσ ¼ 1, with σ min = 0.73, σ max = 1.62, and P(σ) = Aσ −3 , with A a normalization constant 9 .
The dynamical state of the i − th particle is given by its position r i and by the orientation e i of the self propelling force that, in two spatial dimension, is parametrized by the angle θ i , i.e., e i ¼x cos θ i þŷ sin θ i , withx andŷ the unit vectors of the x and y axis, respectively. The overdamped equations of motion for the i − th disk read with 〈η i 〉 = 0 and hη i ðtÞη j ðsÞi ¼ 2 τ δ ij δðt À sÞ. The term ζ i is a thermal noise, i.e., hζ α i i ¼ 0, and hζ α i ðtÞζ β j ðsÞi ¼ 2μk B Tδ ij δ αβ δðt À sÞ. The force F i = ∑ j≠i f ij , with f ij ¼ À∇ r i vðr ij Þ. We perform numerical simulations with k B T = 0.01, μ = 1, v 0 = 1. As control parameter we adopt the persistence length ℓ = τv 0 /〈σ〉 and the packing fraction ϕ = NA s /A, with A s = π〈σ〉 2 . The unit of length is fixed by 〈σ〉 = 1, time is measured in units of τ, and energy in units of ϵ = 1. The system is enclosed in a square box of side L = 60〈σ〉. For exploring the phase diagram, we have performed numerical simulations for N ∈ [15 2  Instantaneous and inherent configurations. The inherent structures have been obtained minimizing the mechanical energy Φ ¼ 1 2 ∑ i≠j vðr ij Þ. Energy minimization is performed using the FIRE algorithm 73 . For gaining insight into the stability of inherent and instantaneous configurations, we have computed the normal modes by evaluating the 2N eigenvalues of the Hessian matrixλ κ , with κ = 1, .., 2N. The computations have been done using Python NumPy linear algebra functions 74 . The density of states of instantaneous and inherent configurations have been computed in a standard way, by introducing the eigenfrequency ω κ ¼ ± ffiffiffiffiffiffiffi jλ κ j p , with negative frequency corresponding to negative eigenvalues (that populate the spectrum only in the case of instantaneous configurations). The density of states reads DðωÞ ¼ N À1 ∑ κ δðω À ω κ Þ, with N ¼ 2N À 2.
Observables. We indicate with hOi s the average of the observable O over independent runs. We indicate with hOi t time-averaging in stationary states. As order parameter describing the global and local properties of the system, we compute the packing fraction field ϕ(x, y, t) and the hexatic field ψ 6 (x, y, t). ϕ(x, y, t) has been obtained discretizing the simulation box into a lattice of linear size ζ = 4〈σ〉. We can thus define a probability distribution function PðϕÞ of the local packing fraction. In the MIPS region homogeneous density profiles become unstable and thus PðϕÞ develops a double-peaked structure.
The hexatic field at time t can be computed starting from its microscopic definition and thus performing the Voronoi tesselation of the particle centers and defining the hexatic order parameter ψ i 6 of the particle i as with N i the number of Voronoi neighbors to the cell i. The angle θ ik is individuated by the two cell centers i and k. Again, we compute the distribution Pðjψ 6 jÞ discretizing the simulation box into a lattice of linear size ζ. We obtain additional information on the structural properties of the system measuring the pair distribution function and the g 6 (r) correlation function defined as g 6 ðrÞ ¼ hψ 6 ðr 0 Þψ Ã 6 ðr 0 À rÞi s;t;r 0 : ð7Þ As dynamical observables we measure the self-part of the Intermediate Scattering Function F self (q, t), and the time-correlation function of the hexatic order parameter C ψ (t) 37 where the wave vector q = (q x , q y ) satisfies the periodic boundary conditions, i.e., q x;y ¼ 2π L ðn x ; n y Þ, with n x,y = 0, ±1, ±2, ... (and avoiding the combination n x = n y = 0). The sample-to-sample fluctuations of F self (q, t) provides a measure of dynamical heterogeneity through the four-point dynamical susceptibility χ 4 (q, t). The position of the peak χ 4 (q, t) defines the typical time scale t = τ 4 of dynamical heterogeneity. Map of displacements field has been computed as Δr(x, y, τ 4 ) 40 . We also measure the relaxation time of shape fluctuations using C ψ (t) defined through the correlation function To quantify the similarity between stationary and inherent configurations described through opportune coarse-grained variables x inst and x IS (for instance, x = |ψ 6 |), we have computed the Kullback-Leibler divergence between the probability distributions P inst ðxÞ and P IS ðxÞ that is defined as D KL ½P inst jP IS ¼ Z dx P inst ðxÞ log P inst ðxÞ P IS ðxÞ : ð10Þ

Data availability
The data that support the findings of this study are available from the corresponding author on reasonable request.

Code availability
The code is available from the corresponding author upon reasonable request.