Fundamental differences between glassy dynamics in two and three dimensions

The two-dimensional freezing transition is very different from its three-dimensional counterpart. In contrast, the glass transition is usually assumed to have similar characteristics in two and three dimensions. Using computer simulations, here we show that glassy dynamics in supercooled two- and three-dimensional fluids are fundamentally different. Specifically, transient localization of particles on approaching the glass transition is absent in two dimensions, whereas it is very pronounced in three dimensions. Moreover, the temperature dependence of the relaxation time of orientational correlations is decoupled from that of the translational relaxation time in two dimensions but not in three dimensions. Last, the relationships between the characteristic size of dynamically heterogeneous regions and the relaxation time are very different in two and three dimensions. These results strongly suggest that the glass transition in two dimensions is different than in three dimensions.

I n two-dimensional (2D) solids, thermal fluctuations destroy crystalline order, displacement correlations increase logarithmically and density correlations decay according to power laws 1,2 . However, there can be long-range bond-orientational order in 2D 2 . The transition from the 2D fluid phase to the solid phase can occur in two steps with an intermediate phase characterized by an exponential decay of the density correlations and a power-law decay of the bond-orientational correlations 1,3 . In contrast, in three-dimensional (3D) solids, fluctuations do not destroy crystalline order 4 , and long-range translational and rotational order emerge together at the freezing transition.
Despite these differences between 2D-and 3D-ordered solids, the formation of an amorphous solid on supercooling a fluid, that is, the glass transition, is generally assumed to have similar characteristics in 2D and 3D 5 . This assumption is reflected in the trivial dimensional dependence of most glass transition theories 6 .
We show that structural relaxation of supercooled fluids in two dimensions is different than in three dimensions. While we find the transient localization often associated with glassy dynamics in three dimensions, we do not find any transient localization in two dimensions if we simulated systems large enough to remove any finite size effects. Furthermore, the temperature dependence of the bond-orientational correlation time is decoupled from that of the translational relaxation time in two dimensions, but these relaxation times have very similar temperature dependence in three dimensions. Along with these differences in structural relaxation, we also find that the characteristic size of regions of correlated mobility, dynamic heterogeneities, increases faster with the structural relaxation time in two dimensions than in three dimensions, and these regions are more ramified in two dimensions than three dimensions. Last, we show that the structural relaxation and heterogeneous dynamics depends on the underlying dynamics in two dimensions.

Results
Structural relaxation. To demonstrate the differences between glassy dynamics in two and three dimensions, we focus on two closely related glass-forming fluids: the 3D 80:20 binary Lennard-Jones system introduced by Kob and Andersen 7 , and its 2D variant that has the same interaction potentials but a 65:35 composition to avoid crystallization 8 . To simulate the relaxation in these systems, we used the standard Newtonian dynamics 9 . We also simulated the 2D system using Brownian dynamics 9 and we comment on the differences between these two dynamics. See Methods for the simulation details and the reduced units used to present the results. We examined three other 2D glass formers and one additional 3D glass former (see Methods for details of the systems). We present some results for the additional 2D glass formers in the main text and in Supplementary Material. The figures showing results for the three additional glass formers are labelled with the labels given in Methods, and figures without labels show results for the Kob-Andersen mixtures.
In 3D, the dominant feature in the dynamics of deeply supercooled glass-forming fluids is transient localization of individual particles 6 , which is illustrated in the inset to Fig. 1a. The transient localization results in characteristic plateaus of the self-intermediate scattering function, F s ðk; tÞ ¼ N À 1 h P n e ikÁ½rnðtÞ À rnð0Þ i (|k| ¼ k ¼ 7.2 in 3D and 6.28 in 2D), shown in Fig. 1a, and the mean square displacement, dr 2 ðtÞ h i¼ N À 1 h P n ½r n ðtÞ À r n ð0Þ 2 i, shown in Fig. 1b. The plateaus extend to longer and longer times on approaching the glass transition. Similar plateaus are observed in the collective  scattering function, Fðk; tÞ ¼ N À 1 h P n;m e ikÁ½r n ðtÞ À r m ð0Þ i, which describes relaxation of the density field (not shown), and in the correlation function quantifying bond-orientational correlations in 3D, C Q (t), shown in Fig. 1c (see Methods for the definition of C Q ). Qualitatively similar slowing down of the translational and bond-orientational relaxation in 3D glass-forming fluids is analogous to the simultaneous appearance of translational and rotational long-range order in 3D crystalline solids.
The transient localization observed in 3D glassy dynamics is absent in 2D, as showed in the inset to Fig. 1d. Correspondingly, there is no intermediate time plateau in the self-intermediate scattering function in 2D (Fig. 1d). The final decay of F s (k; t), which in 3D is well described by a stretched exponential, is replaced by a very slow decay in 2D. The intermediate time plateau in the mean square displacement observed in 3D is replaced by an extended sub-diffusive regime in 2D (Fig. 1e).
However, an intermediate time plateau is observed in the correlation function quantifying bond-orientational correlations in 2D, C C (t) ( Fig. 1f; see Methods for the definition of C C ). Shown in Fig. 3a-c are F s (k;T) and C C (t) for three additional 2D glass formers and they behave similarly. Qualitatively different behaviour of the translational and bond-orientational relaxation in 2D glass-forming fluids is analogous to the absence of the translational and the presence of the bond-orientational long-range order in 2D solids 1 .
To quantify decoupling between translational and bondorientational relaxation, we compare the temperature dependence of the relaxation times characterizing F s (k;t) and C (Q,C) (t), where Q and C refer to 3D and 2D correlation functions, respectively. We define the translational relaxation time t a through the relation F s (k;t a ) ¼ e À 1 and the bond-orientational relaxation time t y through C (Q,C) (t y ) ¼ e À 1 . At the highest temperatures,  the ratio t y /t a is less than one for both the 3D and the 2D glass former (Fig. 2a). However, this ratio stays approximately constant with decreasing temperature for the 3D glass former, but grows monotonically for the 2D glass former. In Fig. 3a-c, we show F s (k;t) and C C (t) for three other 2D glass formers, which demonstrates that the decoupling of the temperature dependence of the translational and the bond-orientational relaxation time is general feature of 2D glassy dynamics. In addition, in Fig. 2b,c, we show that the final translational and orientational relaxation satisfies the time-temperature superposition in 3D but not in 2D, and we show corresponding figure for the 32:68 mixture (see Methods) in 3D in Supplementary Fig. 1 and for 2D in Supplementary Fig. 2. Figure 2c clearly demonstrates the decoupling of the temperature dependence of the translational and bond-orientational relaxation times in 2D.
Dynamic heterogeneities. The non-exponential decay of F s (k;t) is frequently attributed to the emergence of domains, referred to as dynamic heterogeneities, in which the relaxation is spatially correlated and significantly different (faster or slower) than the average relaxation. While we find non-exponential decay in F s (k;t) for 3D and 2D glass formers, the nature of the decay is very different and this difference is mirrored by differences in the heterogeneous dynamics. Shown in Fig. 4a-d are displacement maps showing the centre of a four million particle simulation in 2D at T ¼ 0.45. The maps are created by colouring the particles, whose position is shown at t ¼ 0, according to the magnitude of their displacements |r n (t) À r n (0)| between times 0 and t. The red particles have moved a distance equal to or greater than the diameter of a larger particle. There are large domains of particles that have moved less than a particle diameter even at t ¼ 10,000.
Considering the large dynamically heterogeneous regions in Fig. 4a-d, it is unsurprising that we also find large finite size effects. Shown in Fig. 4e is F s (k;t) calculated for different size systems at the same temperature as shown in Fig. 4a-d. A plateau reminiscent of the plateau in 3D systems is present for the smaller systems but gradually disappears with increasing system size. Similar finite size effects are also evident in the mean square displacement (Supplementary Fig. 3) and the inherent structure dynamics ( Supplementary Fig. 4).
To quantify dynamic heterogeneity shown in Fig. 4a-d, we use a four-point structure factor S 4 (q;t) 15 constructed from overlap functions w n (a;t) ¼ Y(a À |r n (t) À r n (0)|), where Y( Á ) is Heaviside's step function. The parameter a is chosen such that N À 1 h P n w n ða; t a Þi % F s ðk; t a Þ, which results in a ¼ 0.25 in 3D and a ¼ 0.22 in 2D. To characterize the slow domains, we calculate S 4 ðq; tÞ ¼ N À 1 h P n P m w n ða; tÞw m ða; tÞe iqÁðrnð0Þ À rmð0ÞÞ i (note that w n (a;t) restricts the sums over the particles that moved less than a over a time t). The characteristic size of dynamically heterogeneous regions is quantified through the dynamic correlation length x 4 (t), which is determined from fitting S 4 (q; t) for small q to the Ornstein-Zernicke form w 4 (t)/{1 þ [qx 4 (t)] 2 }. Here w 4 (t) is the dynamic susceptibility, which characterizes the overall strength of the dynamic heterogeneity.
In Fig. 5a, we show the correlation between the translational relaxation time, t a , and the dynamic correlation length calculated at t a , x 4 (t a ), for the 3D and 2D glass-forming fluids. While for the 3D system we find that a power law is a poor description for an extended range of t a , and a better description is x 4 (t a )B[In(t a /t 0 )] 2/3 (red line in Fig. 5b), we find that a power law x 4 ðt a Þ $ t b a with b ¼ 1.0 ± 0.1 describes the full range of results well for the 2D system. We show results for an additional glass former in 2D and 3D, the 32:68 mixture, in Supplementary  Fig. 5. Note that similar power-law behaviour was observed in simulations of 2D granular fluids 10 . In Fig. 5b, we show that the relationship between the dynamic susceptibility and the dynamic correlation length is fundamentally different in 3D and 2D. For 3D systems, w 4 (t a )Bx 4 (t a ) 3 at low temperatures, which implies compact dynamically heterogeneous regions. For 2D systems, we observe w 4 (t a )Bx 4 (t a ) 1.5 , which suggests more ramified dynamically heterogeneous regions, see Fig. 4d.
Dependence on the microscopic dynamics. Last, we discuss the dependence of the long-time relaxation in 2D on the underlying microscopic dynamics and two important consequences. In 3D, an important finding is that the long-time dynamics does not depend on the microscopic dynamics; the same long-time dynamics has been observed in simulations using Newtonian 7 , stochastic 11 , Brownian 12 and Monte Carlo 13 dynamics. This result can be rationalized within the mode-coupling approach 14 . Surprisingly, we find that in 2D the long-time dynamics is quite different in the case of microscopic Newtonian and Brownian dynamics. The results corresponding to the those shown in   Fig. 6 for the Brownian case. Notably, the decay of F s (k;t) is strikingly different for the Brownian simulations than for the Newtonian simulations. Importantly, the temperature dependence of the translational relaxation time is also decoupled from orientational relaxation time in the case of Brownian dynamics (Fig. 2a) but the ratio t y /t a is not as large for Brownian dynamics than for Newtonian dynamics. In addition, we find a power-law relationship x 4 ðt a Þ $ t b a between the dynamic correlation length and the relaxation time but with b ¼ 0.36±0.05, which is a different exponent than obtained for Newtonian dynamics, see Fig. 5a. However, we find that the relationship between the strength of the dynamic heterogeneity and the dynamic correlation length in 2D is the same for Brownian and Newtonian dynamics, see Fig. 5b. The latter two results show that the universality of the relationships between the relaxation time and properties of heterogeneous dynamics that we found in 3D 15 is absent in 2D. Furthermore, a full description of heterogeneous dynamics in 2D must also include the influence of the microscopic dynamics, and descriptions solely in terms of the structure or the potential energy landscape are not sufficient in 2D.

Discussion
Glassy dynamics in 2D and in 3D are profoundly different. We presented detailed results for one glass former and we verified that the features of the translational relaxation and dynamic heterogeneity are qualitatively the same for three additional 2D glass formers. Our results call for a re-examination of the present glass transition paradigm in 2D. We note that there is currently no theoretical framework that accounts for the different dynamics observed in the 2D glass-forming systems. However, we note that the dynamic picture of the random first-order transition theory breaks down for dimensions less than two, and has been described as marginal for two dimensions 16,17 . Moreover, insights gained from theoretical analysis of the 2D glassy dynamics and glass transition might shed light onto slow dynamics and the glass transition in 3D. It will also be interesting to investigate if the differences between 2D and 3D glassy dynamics are observable for glass-forming fluids in confinement and at interfaces or surfaces, that is, for quasi 2D systems.

Methods
Simulations. We simulated binary mixtures of Lennard-Jones particle in two and three dimensions. The interaction potential is V ab (r) ¼ 4E ab [(s ab /r) 12 À (s ab /r) 6 ] where E BB ¼ 0.5E AA , e AB ¼ 1.5E AA , s BB ¼ 0.88s AA and s AB ¼ 0.88s AA . The results are presented in reduced units where s AA s is the unit of length and E AA the unit of energy. The unit of time for the Newtonian dynamics simulations is ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi s 2 m=E AA p and the mass m is the same for both species. The Newtonian dynamics simulations 9 were performed using LAMMPS 18 for the 2D and 3D simulations and HOOMD blue for the 2D simulations 19 . The LAMMPS simulations were run in an  ARTICLE NVE ensemble but there is significant energy drift for the HOOMD blue simulations for the lowest temperatures. Therefore, we ran the HOOMD blue simulations using an NVT Nosé-Hoover thermostat with a coupling constant t ¼ 10. We ran at least one LAMMPS NVE simulation at every temperature to make sure that the conclusions did not depend on the thermostat. All the results are averages over four or more production runs. The equations of motion for the Brownian dynamics simulations 9 are r n (t) ¼ g À 1 F n (t) þ g n (t), where g ¼ 1 is the friction coefficient, F n (t) is the force on particle n at time t and g n is a random noise term. The random noise satisfies the fluctuation dissipation relation Z n ðtÞZ m ðt 0 Þ ¼ 2k B Tg À 1 dðt À t 0 Þd nm 1, where 1 is the unit tensor. The unit of time for the Brownian dynamics simulation is s 2 g/E AA . The Brownian dynamics simulations were run using a modified version of LAMMPS and our in house developed code. Inherent structure trajectories were created by quenching the system to the local potential energy minimum using the FIRE algorithm 20 .
We simulated 2D systems of 10,000 particles for TZ0.9 and 250,000 particles for 0.5rTr0.8. At T ¼ 0.45, we studied 4 million particles for the Newtonian dynamics simulations, but 250,000 particles for the Brownian dynamics simulations. We simulated 27,000 particles in 3D using Newtonian dynamics. To check that the results are independent of the system size for each state point for the 2D Newtonian dynamic simulations, we ran 100,489 particle simulations and checked to see if the results agreed with the 250,000 particle simulations. At T ¼ 0.45 they did not agree, and we increased the system size until we found agreement between the 4 million particle system and an 8 million particle system. For the 2D Brownian dynamics simulations, we found agreement between 10,000 particle simulations and 250,000 particle simulations for TZ0.45. For T ¼ 0.4, we found that a 100,489 particle simulation agreed with a 250,000 particle system.
We also examined the translational dynamics, bond-orientational and dynamic heterogeneities for three additional systems in 2D and one additional glass-forming system in 3D. The first system is the one studied in ref. 21 and consists of a 31.67:68.33 binary mixture with the potential V ab (r) ¼ E(s ab /r) 12 . The size ratios are s AB ¼ 1.1s AA and s BB ¼ 1.4s AA . We simulated this system using 250,000 particles in 2D and 100,000 particles 3D. The number density r ¼ 0:719s À 2 AA in 2D and r ¼ 0:719s À 3 AA in 3D. The units for energy is E, length is s AA , time is ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi s 2 AA m=E p , and temperature is E/k B . We denote this system with the label 32:68. The second system was introduced in ref. 22, and consists of an 50:50 mixture of repulsive particles where the potential V ab (r) ¼ E(s ab /r) 12 . The size ratios are given by s AB ¼ 1.2s AA and s BB ¼ 1.4s AA and the number density r ¼ N=L 2 ¼ 0:74718s À 2 AA . We simulated 250,000 particles for this second additional system. The units for energy is E, length is s AA , time is ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi s 2 AA m=E p and temperature is E/k B . We denote this system with the label 50:50.
The third system is the one introduced in ref. 23, which consists of an 50:50 mixture of harmonic spheres with the interaction potential V ab (r) ¼ 0.5E(1 À r/s ab ) 2 for r s ab and V ab (r) ¼ 0 otherwise. The size ratios are given by s AB ¼ 1.2s AA and s BB ¼ 1.4s AA and r ¼ 0:699s À 2 AA . We simulated 250,000 particles for this third additional system. The units for energy is E, length is s AA and time is ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi s 2 AA m=E p . The unit for temperature is 10 À 3 E/k B . We denote this system with the label Harm.
Bond-orientational correlation functions. To measure bond-orientational relaxation times in 2D, we first define C n 6 ðtÞ ¼ ðN n b Þ À 1 P m e i6ynm ðtÞ , where y nm (t) is the angle between particle n and particle m at a time t, N n b is the number of neighbours of particle n, and the sum is over the neighbours of particle n at the time t. The neighbours are determined through Voronoi tessellation 3 . The time dependence of the bond angle correlations was monitored by calculating C C ðtÞ ¼ h P n C n 6 ðtÞ½C n 6 ð0Þ Ã i=h P n j C n 6 ð0Þj 2 i, where * denotes the complex conjugate.
To measure bond-orientational relaxation in 3D, we define Q i lm ðtÞ ¼ ðN i b Þ À 1 P j q lm ½y ij ðtÞ; f ij ðtÞ, where q lm (y,f) are the spherical harmonics 24 and the sum is over the neighbours of a particle i at a time t determined through Voronoi tessellation. Next, we define the correlation function Q l ðtÞ ¼ h½ð4pÞ=ð2l þ 1Þ P i P l m¼ À l Q i lm ðtÞ½Q i lm ð0Þ Ã i. We calculated C Q (t) ¼ Q 6 (t)/Q 6 (0) to monitor the decay of orientational correlations.
We note that the conclusions remain unchanged if we define neighbours as being less than a distance equal to the first minimum of the pair correlation function rather than through Voronoi tessellation.