Results: Unpinned Case

The system studied is composed of N=10,000 polydisperse hard disks with disk-size polydispersity Δ=11% (Materials and Methods). We start by considering the case without an external pinning field, c=0. We focus on the following number densities from ρ=0.92 to ρ=0.97, which correspond to area fractions ranging from ϕ=0.73 to ϕ=0.77. As described in SI Appendix, we extract the correlation length for bond-orientationally ordered regions (SI Appendix, Fig. S1) by fitting the exponential decay of the peaks of the correlation function g6(r)/g(r) (SI Appendix, Fig. S2). We use an exponential function instead of a 2D Ornstein–Zernike function to avoid a priori assumptions on the origin of the growth of the correlation length. The correlation length ξ6 is plotted in Fig. 1D, together with other length scales that we will derive later. The two-body correlation function ξ2 is obtained by fitting with an exponential law the decay of g(r)−1 (SI Appendix, Fig. S3). The results of this fit are also summarized in Fig. 1D. These results confirm that, for polydisperse glass-forming systems, the growth of many-body correlations associated with bond-orientational order is much faster than the growth of two-body correlations (11, 13, 32).

As was shown in refs. 11 and 13, there is a link between dynamic heterogeneities and regions of high hexatic order in the fluid. We provide here further support to this scenario. We measure dynamical length scales by performing event-driven molecular dynamics simulations and extracting the four-point density correlator, with the procedure described in SI Appendix, Figs. S4 and S5. In Fig. 1D, we plot the dynamical correlation length ξ4, showing that indeed ξ4 scales with density like ξ6. To show that there is a causal relation between slow dynamical regions and regions of high hexatic order, in Fig. 1 E and F, we plot snapshots of configurations prepared at ρ=0.92 and ρ=0.97, respectively. For each snapshot we run event-driven molecular dynamics simulations in the isoconfigurational ensemble (Materials and Methods), where 200 trajectories are started from the same initial configuration but with different initial velocities. The degree of structural order is investigated by taking the average hexatic field over these N=200 trajectories (also called isoconfigurational average) at a time t=τα/10, where τα is the structural relaxation time measured through the intermediate scattering function (SI Appendix, Fig. S6). The dynamics is instead investigated through the isoconfigurational average of the displacement |Δri| between t=0 and t=τα, which is also approximately the time at which the heterogeneities are maximum (as measured by the four-point susceptibility shown in SI Appendix, Fig. S4). All disks in the snapshots of Fig. 1 E and F are then grouped into sets of high and low mobility/order, depending on whether their mobility/order is higher or lower than the 50th percentile. We can then identify four different sets of particles: low mobility and high order (white); high mobility and low order (black); low mobility and low order (orange); and high mobility and high order (magenta). Our results show that, for all densities considered in this work, 66% of particles are in the first two sets (33% in each), demonstrating a high degree of correlation between structural ordered regions and immobile regions (or, vice versa, between disordered regions and mobile regions) even at a particle level. Moreover, the remaining two sets (each accounting for the 17% of particles) are located at the interface between mobile and immobile extended regions. Magenta disks are located on the surface and in between black clusters, whereas orange disks are located on the surface and in between white clusters. In other words, disks that are next to a low/high mobility region will also have low/high mobility. Here, we note that the embedded fractal nature of order parameter fluctuations is characteristic of critical fluctuations. In SI Appendix, we also note that a higher degree of correlation can be obtained if the relative displacement |ΔRi| is used instead of the absolute displacement |Δri|. The relative displacement |ΔRi| is defined as the displacement between time t=0 and t=τα, of particle i with respect to its M neighbors, Ri=ri−(1/M)∑jMrj. The coherent motion of a particle with its neighbors does not contribute to the stress relaxation as no bonds are broken in this process. In this case, the first two sets account for 76% of the particles (SI Appendix, Fig. S7).

Results: Pinning

Having characterized the static and dynamic properties of the unperturbed system, we now introduce the pinning field. First, simulations are fully equilibrated with cluster-moves algorithms (Materials and Methods). A representation of the pinning fields is given in Fig. 1 A–C. In the “random-pinning” geometry, Np particles are chosen randomly and pinned, i.e., their position is kept fixed during the course of the simulations. For each density, we introduce pinning fields with concentrations c=0.01, 0.06, 0.10, 0.15, and 0.20, and for each concentration we average over nine different realizations of the fields. In this scheme, the distance between pinned particles is defined only as an average over a broad distribution, as both clusters of pinned particles and extended regions without pinned particles are likely produced. For this reason, random pinning is expected to be more sensitive to finite size effects, as was observed in ref. 19, where it is noted that random pinning can smear out the Kauzmann transition in very small systems. The pinning geometry can also have strong effects on the dynamics (21, 30, 33, 34). To limit the fluctuations in the distance between pinned particles, we also adopt a “uniform-pinning” geometry, where a simple cubic lattice is overlaid to the equilibrated configuration, and the closest particle to each lattice point is pinned. Particles are pinned at the following average distances: a=2.5, 2.75, 3, 4.25, 6, 8, and 10, and each distance is averaged over seven realizations of the field. Finally, in the “cavity-pinning” geometry, all particles outside a cavity of radius R are pinned. Because the static length scales currently accessible to simulations are expected to be small, well within 10σ, simulations with cavity pinning involve a small number of particles, thus requiring extensive average over different realizations of the field (here 100 simulations for each cavity diameter).

Because the pinning field is applied to equilibrium configurations, the static properties should be unchanged with respect to the c=0 case. We check this by computing both positional and hexatic order for different concentrations c. All results are consistent with the c=0 case and the SD between simulations at different c is represented with the error bars for ξ6 in Fig. 1D. Correlation lengths are extracted from all pinning geometries, following the procedure outlined in SI Appendix, Figs. S8–S12. In all cases, the physical idea is to detect the characteristic length (average distance between pinned particles in the random- and uniform-pinning geometries, or the size of the cavity in the cavity geometry), which produces a high localization of the mobile particles, as measured by overlap functions. We plot the random-pinning correlation length ξK,random, the uniform-pinning correlation length ξK,uniform, and the PTS length scale from cavity pinning ξPTS in Fig. 1D (see SI Appendix, Figs. S8–S12, on the details of its estimation). We see that the growth of the PTS length scale, irrespective of the pinning strategy, is significantly slower than the growth of bond-orientational order ξ6, while being comparable to the growth of pair correlations ξ2 in the system. We also confirm that estimating the PTS length from uniform and random pinning through the relation ξPTS(T)∼ξK2(T) still produces a much weaker growth than that of ξ6. The Inset of Fig. 1D shows the scaling between the ξ6 length scale and the relaxation time τα, τα=τ0⁡exp(Dξ/ξ0), where D is a measure of the fragility of the system (10). This scaling also supports a direct connection between the growth of structural correlation and slow dynamics, but its origin is still unclear. A justification of this relation in the context of a theory of topological ordered clusters has been given in ref. 26.

All static length scales can also be obtained by considering coarse-grained variables, dividing the simulation box into smaller boxes of side length l=0.3σ, ensuring that each box can be occupied at most by one disk at any time during the simulation. The problem is then mapped on a set of discrete variables, defined as niα=1 when the center of a disk is in box i in configuration α, and niα=0 otherwise. The average occupancy for a particular realization of the pinning field is defined as ni¯, where the overline denotes an average over thermal fluctuations for a fixed realization of the pinning field. Fig. 2 shows the field ni¯ for three different concentrations of the pinned particles for the case of uniform pinning at ρ=0.92 (top row) and ρ=0.97 (bottom row). Going from the low to high concentration (Left to Right in the figure), the amount of localization progressively increases. Localization occurs when the occupancy field becomes strongly peaked in discrete locations in space, surrounded by regions of negligible probability occupancy. For ρ=0.92 (top row in Fig. 2), the displayed pinning distances a are always bigger than the pinning correlation length (Fig. 1D), and the occupancy field ni¯ is homogeneous. Instead, for ρ=0.97 (bottom row in Fig. 2), regions of high localization progressively appear and form a connected network close to a=3, which is the pinning correlation length at this density (ξK,uniform∼a∼3). So the pinning correlation length ξK is the length scale at which localized regions are close to percolation. In SI Appendix, Fig. S13, we show that, at the level of undercooling that we can reach, the transition from a homogenous state to a localized state is continuous.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Development of ni¯ field with a decrease in a for uniform pinning. The particle density is ρ=0.92 (top row) and ρ=0.97 (bottom row), and the average distance between particles a=10 (Left), a=4.25 (Middle), and a=3 (Right). The color code is mapped to the occupancy probability according to the color bar. The Inset inside some of the panels shows the magnification of the Bottom Left corner of the image of an area of one-eighth of the total image. Distances are in units of the average diameter 〈σ〉.

A visual inspection of Fig. 2 (bottom row) shows that localized particles have a high degree of hexatic order (SI Appendix, Fig. S13), already suggesting that ξK is smaller than the hexatic correlation length ξ6. From Fig. 1D, we know that the pinning length scale ξK is close to the pair correlation length ξ2. Each pinned particle generates an oscillatory perturbation of the ni¯ field, which originates from the two-body static correlations between the pinned particle and the mobile particles in the liquid, and whose extent is thus given by ξ2. The length scales measured in Fig. 1D thus indicate that the localization occurs when the average distance between pinned particles is comparable to the range of two-body correlations, ξ2, below which the number of particle arrangements drastically decreases, and the configurational entropy vanishes. The extent of the regions with high localization of particles is exactly what is being measured by the PTS correlation function.

The results thus show that the growth of the pinning correlation length is similar to the growth of two-body correlations. On the other hand, the growth of bond-orientational order is much faster and clearly decoupled from the pinning correlation length. We also checked that the same is true for coarse-grained quantities. This strongly indicates that the PTS correlation length is not order agnostic but targets the growth of a particular order in the system, that is, the size of the regions where particles are localized due to the pinning field. The growth of these regions follows the growth of two-body correlations: particles are localized due to the perturbation that pinned particles introduce to the ni¯ field, and the length scale of this perturbation is given by two-body correlations. In other words, at least in the density range considered here, the localization transition due to point pinning requires that the average distance between pinned particles is smaller than the two-body correlation length.

Discussion and Conclusions

In this article, we have extracted several static length scales from systems of polydisperse hard disks with polydispersity Δ=0.11, in the range ρ∈[0.92;0.97]. The results confirmed that the length scale associated with bond-orientational order grows more rapidly than the length of pair correlations (11, 13, 32). The use of pinning fields enabled the calculation of the PTS correlation length, showing that it grows only moderately with increasing supercooling, a result that is in agreement with measures of the PTS length in binary mixtures (23). For polydisperse systems, the PTS correlation length is not coupled to that of bond-orientational order, which is directly linked to the dynamical correlation length (10, 11): the growth of the former is considerably slower than the latter. For different glass-forming systems, this suggests that also other forms of order originating from many-body interactions could go undetected by PTS measures.

The PTS length captures a localization transition that occurs in presence of pinned particles. This localization transition originates when the occupancy field, n¯, has extended regions of high probability due to neighboring pinned particles. The perturbation that a single pinned particle produces in the n¯ field is due to “pair correlations.” In absence of strong nonlinear effects, the length scale of the localized regions extends no further than two-body correlations. This is the case in the density interval accessible to our simulations, where pinned particles need to be placed closer than the pair correlation length to produce localized regions in the fluid. A second requirement that our results suggest is that pinned particles should be in positions compatible with high local hexatic order.

The localization transition that occurs with increasing concentration of pinned particles happens continuously in the density range we could access in equilibrium. The results do not rule out the possibility that strong nonlinearities will produce a localization transition that extends beyond pair correlations for higher (but yet unreachable) densities.

To summarize, the PTS length measured by particle pinning simply reflects pair correlation and fails in detecting the correlation of bond-orientational order (more precisely, hexatic order), which intrinsically originates from many-body interactions. Although the PTS length is decoupled from the dynamical correlation length, the hexatic order correlation is strongly coupled to it. This implies that slow dynamics in our system is controlled by the development of hexatic ordering, and not by translational order detected by the PTS correlation. Although the generality of this conclusion needs to be checked carefully, our study suggests that the PTS length is not order agnostic and the growth of the PTS length is not responsible for glassy slow dynamics, at least for our system.

Differently from the fluid-hexatic transition, the results in refs. 11 and 13 have provided evidence that the glass transition in polydisperse hard disks involves a power-law growth of the hexatic correlations that follows the Ising universality class, but, differently from ordinary critical phenomena, is accompanied by a strong divergence [log(τα)∼ξd/2] of the relaxation time. This behavior could originate from frustration and disorder effects, which cause the energy barriers for viscous flow to increase with decreasing temperature (10). Understanding the relation between power-law growth of the ξ6 and the divergence of dynamics represents the main challenge to be addressed in future work.