Transport signatures of fragile-glass dynamics in the melting of the two-dimensional vortex lattice

In two-dimensional (2D) systems, the melting from a solid to an isotropic liquid can occur via an intermediate phase that retains orientational order. However, in 2D superconducting vortex lattices, the effect of orientational correlations on transport, and their interplay with disorder remain open questions. Here we study a 2D weakly pinned vortex system in amorphous MoGe films over an extensive range of temperatures ($\bm{T}$) and perpendicular magnetic fields ($\bm{H}$) using linear and nonlinear transport measurements. We find that, at low fields, the resistivity obeys the Vogel-Fulcher-Tamman (VFT) form, $\bm{\rho(T)\propto\exp[-{W}(H)/(T-T_0(H))]}$, characteristic of fragile glasses. As $\bm{H}$ increases, $\bm{T_0(H)}$ is suppressed to zero, and a standard vortex liquid behavior consistent with a $\bm{T=0}$ superconducting transition is observed. Our findings, supported also by simulations, suggest that the presence of orientational correlations gives rise to a heterogeneous dynamics responsible for the VFT behavior. The effects of quenched disorder become dominant at high $\bm{H}$, where a crossover to a strong-glass behavior is observed. This is a new insight into the dynamics of melting in 2D systems with competing orders.

In two-dimensional (2D) systems, the melting from a solid to an isotropic liquid can occur via an intermediate phase that retains orientational order. However, in 2D superconducting vortex lattices, the effect of orientational correlations on transport, and their interplay with disorder remain open questions. Here we study a 2D weakly pinned vortex system in amorphous MoGe films over an extensive range of temperatures (T ) and perpendicular magnetic fields (H) using linear and nonlinear transport measurements. We find that, at low fields, the resistivity obeys the Vogel-Fulcher-Tamman (VFT) form, ρ(T ) ∝ exp [ − W (H)/(T − T 0 (H))], characteristic of fragile glasses. As H increases, T 0 (H) is suppressed to zero, and a standard vortex liquid behavior consistent with a T = 0 superconducting transition is observed. Our findings, supported also by simulations, suggest that the presence of orientational correlations gives rise to a heterogeneous dynamics responsible for the VFT behavior. The effects of quenched disorder become dominant at high H, where a crossover to a strong-glass behavior is observed. This is a new insight into the dynamics of melting in 2D systems with competing orders.
Thermal melting of 2D crystalline solids has been investigated in a variety of systems, including colloids, electrons on the surface of liquid helium, raregas atoms on substrates such as graphite, liquid crystal films, dust plasmas, and vortices in thin superconducting (SC) films in a magnetic field [1,2]. The melting transition is generally understood to be driven by the proliferation of topological defects, as described by the Berezinskii-Kosterlitz-Thouless-Halperin-Nelson-Young (BKTHNY) theory [3][4][5]. According to BKTHNY, for weak enough quenched disorder the transition from the 2D (weakly pinned) solid to the liquid phase occurs via an intermediate phase called hexatic. In the hexatic phase, free dislocations appear, breaking the quasi-longrange positional order, but preserving the hexagonal orientational one. By increasing T further, free disclinations form and a fully isotropic liquid is established. This melting sequence has indeed been detected by scanning-tunneling-spectroscopy (STS) imaging of vortices in 2D superconductors [6][7][8][9], but the morphology of dislocations in real systems may be more complicated and depend on microscopic details [1,2]. For example, an STS study of W-based SC films reported a coexistence of liquid with hexatic, as well as with smectic-like (striped) regions of vortex arrangements, and suggested that this might be a feature of melting of 2D solids formed by linearlike units, including vortices in Abrikosov lattices [6]. For weak pinning and intermediate H, numerical studies also support [10] the appearance of grain boundaries, formed by chains of dislocations, which separate the domains of topologically ordered vortex lattice. The presence of competing orders may give rise to metastable states and the associated slow dynamics in many systems [11]. An STS study on amorphous MoGe (a-MoGe) thin films did observe [9] a strong suppression of the vortex diffusivity in the presence of hexatic correlations compared to the isotropic liquid, but the nature of the dynamics near the melting transition and its effect on electrical transport have not been explored. Another open question is the evolution of transport properties with H in the presence of orientational correlations. The random potential generated by quenched impurities is indeed coupled to the vortex density, so that the increase of H increases the effective disorder [12].
Here, we perform extensive magnetotransport measurements on a-MoGe thin films similar to those used in the STS studies that revealed the presence of hexatic correlations at low H [8,9]. Our key experimental finding is the vanishing of the resistivity, in the VFT fashion, at a nonzero temperature T 0 (H), for fixed low fields H < 9 T. The VFT law [13][14][15] generally describes the behavior of the socalled fragile liquids above the glass transition, and it is usually attributed to the emergence of dynamical heterogeneities [1]. As H increases, we find that T 0 (H) is suppressed to zero and remains zero over a finite range of H (9 < H 12.5 T). In this highfield range, both the observed Arrhenius dependence (T 0 = 0) of the linear resistivity and nonlinear transport (i.e. voltage-current characteristics V -I) are consistent with the disorder-induced freezing of an isotropic vortex liquid to an amorphous vortex glass phase at T = 0. Our finding of the fragile-glass behavior at low H, i.e. in the regime where disorder is not dominant, strongly suggests that the slow vor-tex dynamics described by the VFT law results from dynamical heterogeneities induced by emergent orientational correlations as the melting transition is approached. To elucidate this connection further, we complement our experimental work with Monte Carlo simulations of the 2D XY model in the presence of a transverse magnetic field, analyzing both the static and the dynamic properties of the system. We find that the orientational order leads to a caging effect that suppresses vortex diffusivity, analogous to caging effects induced by disorder in fragile glass-forming liquids that exhibit VFT behavior [1]. Indeed, we identify several signatures of the heterogeneous nature of the vortex dynamics in the presence of hexatic correlations. Our study of transport thus reveals a fragile-glass-like dynamics associated with the thermal melting of a weakly pinned 2D vortex lattice, and a crossover to a strong-glass (T 0 = 0) behavior at higher H, resulting from the interplay of orientational correlations and disorder. The data are shown for sample 1. a, ρ vs H for several T , ranging from 0.017 K to 14 K. b, ρ(T ) at selected H extracted from the data in a. The dashed lines guide the eye. c, ρ(T ) fitted with the Arrhenius (red dashed line) and VFT (blue dotted line) forms for H = 5 T; the fits are performed for the data points within the two red arrows. ρ at the lowest T drops much faster than predicted by the Arrhenius function, but such a fast drop is captured well by the VFT fit; herẽ T0 = (3.42 ± 0.01) K [the fit to Eq. (3) gives T0 = (3.37 ± 0.01) K]. In fact, the VFT fit describes the data even at somewhat higher T , beyond the red arrow. d, For H = 11 T, all the low-T data are described by the Arrhenius form (or VFT form with T0 ≈ 0 K), as shown by the red dashed curve.

Linear resistivity measurements
Our samples are 22 nm thick a-MoGe films with SC transition temperatures T c ∼ 7.6 K at zero field (see Methods), and low normal-state sheet resistances R s ≈ 78 Ω indicative of a very weak disorder [9,17].
Here T c is defined as the temperature at which the linear resistance R ≡ lim I→0 V /I (i.e. ρ) becomes zero, that is, falls below the experimental noise floor (see also Methods). The characteristic length scales for vortex distortions parallel to the applied H and caused by thermal fluctuations or pinning are of the order of several µm [18,19], i.e. much larger than the sample thickness. Therefore, the vortex lattice (VL) in these films is indeed 2D, although the SC state is 3D (thickness > ξ, where ξ ∼ 5 nm is the SC coherence length [18]).
Measurements of ρ(H) at fixed T (Fig. 1a) were used to determine T c (H) (i.e. the corresponding field H c for a given T ) and the upper critical field H c2 (T ), which we define as the magnetic field where ρ = 0.95 ρ N at a given T (ρ N = 0.17 mΩ cm is the normal-state resistivity). We note that, since T c (H) depends on the experimental resolution, it does not necessarily allow one to determine the real boundary between the liquid and solid phases, nor the possible transition to an orientational liquid phase.
To explore the melting of the VL, we extract ρ(T ) curves at fixed H ( Fig. 1b; also Fig. 4a). A rapid, orders-of-magnitude change of ρ with T , observed for fields below about 12 T and with no sign of saturation at low T , suggests an exponential ρ(T ) dependence. Indeed, the Arrhenius law, is often used to describe the linear resistivity in the vortex liquid regime at low enough T [20]. In this, thermally assisted flux flow (TAFF) picture, vortices move collectively and overcome pinning barriers via thermal excitations when T U (H). At the same time, although ρ is nonzero, V -I remains nonohmic for I = 0 [12]. We note that Eq. (1) assumes that ρ is finite at all nonzero T . However, when the superfluid phase sets in, which corresponds to a solid phase for the VL, ρ vanishes, so deviations from the behavior (1) should be observed whenever the critical temperature is finite and the transition is not strongly first order.
We have performed fits of the ρ(T ) data, for fixed H, according to Eq. (1). Two examples are shown in Fig. 1, one for low fields (H = 5 T, panel c) and another for high fields (H = 11 T, panel d). While at high H the Arrhenius fit works extremely well down to T c , at low fields we systematically find significant deviations (see also Fig. 5), with a faster suppression of ρ than predicted by Eq. (1). Interestingly, such deviation can be captured very well by the VFT law, HereW (H) is independent of T , andT 0 (H) provides the temperature where ρ is expected to vanish. Since we are interested in the vortex diffusivity, we have also performed fits to ρ( is the quantum of magnetic flux attached to a single vortex, h is Planck's constant, e is electron charge), k B is the Boltzmann constant, and the vortex diffusion coefficient For both (2) and (3), the fitting parameters were extracted using global minimization, and typically found to be the same within error (see, e.g., Fig. 1c and Fig. 5), that is,Z and Z are the same within the proportionality constant between ρ and D v ; this is consistent with the exponential term dominating ρ(T ). Hence, hereafter we present only the VFT fitting parameters obtained using Eq. (3).
The results for T 0 (H) are shown in Fig. 2, along with the values of T c (H) and H c2 (T ) (see also Fig. 4b; Fig. 6 shows all fitting parameters). We note that, when comparing the Arrhenius and the VFT fits, we paid attention not only to the quality of the fit, but also to the range of temperatures where either one was effective; those ranges are also shown in Fig. 2. As a result of such analysis, we find that for fields below H * 9 T the VFT fit gives a significantly better description of the data. With increasing H, however, T 0 decreases gradually, and vanishes near H * 9 T. For 9 H 12.5 T, ρ(T ) curves are well described by the Arrhenius fits, consistent with T 0 = 0 (Fig. 7). Furthermore, in this regime, we find that U (H) = U 0 ln(H 0 /H) [ Fig. 8; U 0 = (31.5 ± 0.4) K and H 0 = (12.4 ± 0.5) T], as expected from the TAFF model and the logarithmic vortex-vortex interactions in 2D [20].

Phase diagram
The behavior observed in the high-field (H > H * , Fig. 2) is, therefore, consistent with the thermally-activated collective motion of vortices in the presence of strong pinning, similar to previous studies of disordered SC films, including a-MoGe [21,22]. This conclusion is further supported by our measurements of the differential resistance (dV /dI) versus dc current bias I dc at fixed H (see Methods). Indeed, here the V -I characteristic at low enough T remains non-ohmic for I dc = 0 and dV /dI increases with I dc (Fig. 9a, b). Such nonlinear transport is expected from the motion of vortices in the presence of disorder, i.e. it T0(H) (blue dots) are obtained from the VFT fits; the error bars reflect the standard errors of the fits. In the shaded blue "VFT" region, ρ(T ) obeys the VFT law for a fixed H 9 T. For H 9 T, the resistivity above Tc, in the shaded orange "Arrhenius" region, can be only fitted with the Arrhenius law (1). The purple dot-dashed line shows the high-T extent of the VFT and Arrhenius fits. Green triangles: Hc2(T ), i.e. the upper critical field; here ρ reaches 95% of its normal-state value. Red diamonds: Tc(H), where ρ drops below the experimental noise floor. Error bars in Tc(H) and Hc2 represent the uncertainty in defining the magnetic field corresponding to each quantity within our experimental resolution.
is a signature of a viscous vortex liquid [12]. As T increases, nonlinear behavior is no longer observed (Fig. 9c, d). At T < T c , where ρ drops below our noise floor, a finite I dc (i.e. a critical current) is then needed to depin a measurable number of vortices (Fig. 9e). Since ρ(T ) in the liquid phase can only be fitted with the Arrhenius law (1), our transport results suggest that in the high-field regime an isotropic vortex liquid freezes into an amorphous vortex glass as T → 0, i.e. that the solid phase is only realized at T = 0.
In the normal state (H > H c2 ), V -I characteristics are ohmic (Fig. 9b, e), as expected. In the low-field (H < H * , T 0 = 0, T > T c ) regime ("VFT" in Fig. 2) we also find Ohmic behavior (Fig. 9f). However, since this regime is measured only at relatively high T , in analogy with the Arrhenius region we speculate that any nonlinear transport that might be present for small I dc at lower T is experimentally inaccessible. Nevertheless, this leads us to the question about the origin of the observed VFT behavior.
Within the context of glasses, one can define a temperature T 0 , where the relaxation time diverges and the configurational entropy vanishes [1]. T 0 is below the temperature T g where the dynamical glass transition occurs, the value of which depends on the sensitivity of the probe. The VFT fit allows one to overcome such an experimental upper bound and to infer T 0 , that is, the temperature scale at which the residual entropy is comparable to that of the ordered state. Within the context of our experiment, the role of T g is then played by T c , where the resistance drops below the measurable threshold, while T 0 , where the diffusion coefficient (or ρ) vanishes, corresponds to the temperature where the transition to the truly superfluid phase occurs. In that case, we expect that the "VFT" region above T c shows the extent of the dynamically heterogeneous vortex liquid in the (T, H) phase diagram (Fig. 2). Interestingly, analogies to glasses were used to propose the VFT law (2) to describe the slowing down and freezing of the strongly disordered 3D vortex matter [23]. However, in our a-MoGe films the disorder-dominated high-field regime (T 0 = 0) is described, in contrast, by the TAFF model (1). Thus our results strongly suggest that disorder is not the main origin of the VFT suppression of the vortex diffusivity observed at lower fields. To explore the possibility that the presence of orientational correlations may give rise to dynamical heterogeneities and the VFT behavior, we performed the following simulations.

Monte Carlo simulations
We performed Monte Carlo (MC) simulations on the 2D XY model in a transverse field. Its Hamiltonian reads: where θ i represents the SC phase of the condensate, J the effective Josephson-like interaction between nearest neighboring sites, and F µ i is the Peierls phase resulting from the minimal substitution prescription, The intensity of the magnetic field Bẑ = ∇ × A can be expressed in terms of the quantum-flux fraction f passing through a unitary plaquette f = Ba 2 /Φ 0 , where a = 1 is the lattice spacing.
Within the model (4), vortices appear as topological excitations of the phase variable θ i , allowing for the direct characterization [24][25][26][27][28][29] of the static properties of the solid to liquid transition, via computation of the phase rigidity. On the other hand, dynamical effects have been mainly studied via effective models where vortices are mapped into individual particles [10,[30][31][32]. Here we show how model (4) allows us to address both aspects. Indeed, the solid phase can be identified via the superfluid response, and the dynamical properties of vortices can be characterized by tracking the diffusion of each individual vortex in time at a given temperature.
To establish the superfluid phase, we computed the superfluid stiffness J s , namely the global phase rigidity of the condensate (see Methods). At the same time, to establish the orientational order of the VL, we computed the six-fold orientational order parameter G 6 , where the sum runs over the N v vortices of the lattice, ψ 6j indicates the local orientational order relative to the j−th vortex (see Methods for details), and · · · denotes the thermal average and the average over 10 independent numerical simulations. The temperature dependence of J s and G 6 is shown in Fig. 3a, along with prototypical images of the corre-sponding VL; the temperature is expressed in units of J/k B . We can identify three distinct phases: a low-T SC phase, where the VL is a pinned solid with a complete hexagonal order (G 6 1); a disordered non-SC phase at high T , where the VL has melted into an isotropic liquid (G 6 0); and an intermediate phase which is a liquid (J s = 0) but with a persistent orientational order of the VL (G 6 = 0) (see also Supplementary information and Fig. 10). The isotropic to hexatic liquid critical temperature T hex in Fig. 3a has been identified from the orientational susceptibility χ 6 (see Methods and Fig. 11). Due to the small size of our system we cannot discriminate between a floating solid [27,31] and a hexatic liquid based on the static properties. Nonetheless, the dynamic features of the vortices show all the typical fingerprints of the hexatic phase, as they have been identified in 2D soft-colloidal systems [33][34][35][36].
To study vortex dynamics we tracked in time each individual vortex, 1. Above T0, the quasi long-range positional order is lost due to the appearance of dislocations, formed by bound pairs of disclinations with opposite sign (marked with blue and orange dots in the VL snapshots). As T further increases above T hex , isolated disclinations appear, leading to a fully isotropic liquid with vanishing orientational order. Notice that, due to the small lattice size, G6 remains finite and disclination pairs always appear nearby. and we computed the vortex mean-square displacement ∆r 2 (t) at several temperatures, as shown in the inset of Fig. 3b (also Supplementary information and Fig. 9).
At high T , ∆r 2 (t) directly crosses over from the short-time subdiffusive regime (due to the presence of the numerical square grid) to the typical long-time diffusive behavior, where ∆r 2 (t) ∼ D v t, with D v the vortex diffusion constant. However, at the verge of the isotropic-to hexatic-liquid transition, an additional subdiffusive regime appears with a strong suppression of ∆r 2 (t) at intermediate time-scales that is the hallmark of dynamical heterogeneities; here this is due to the caging mechanism provided by the finite orientational order in the hexatic phase (see also Supplementary Video S1). The signatures of the dynamical heterogeneities persist at longer time scales, with a marked reduction at these temperatures of the asymptotic diffusion coefficient, defined as D v = 1 4 lim t→∞ ∆r 2 (t) /t. The resulting temperature dependence of D v (T ) is shown in Fig. 3b. It is clear that the Arrhenius law (continuous red line) strongly deviates from D v (T ) at low T , in contrast to the VFT law (dashed blue line), which provides a good description of D v (T ), similar to our experimental findings in the lower-field regime.
The most remarkable observation is that the temperature T 0 extracted from the VFT fit in Fig. 3b almost coincides with the temperature where J s vanishes (Fig. 3a), showing an internal consistency between the information available from the vortex dynamics and the static properties of the system in the description of the transition from solid to liquid.

Discussion
Our study has revealed similarities between the thermal melting of a 2D VL in weakly pinned a-MoGe films and the behavior of fragile glass-forming liquids. The results of our numerical work strongly support the existence of a glassy dynamics in the absence of strong disorder. Using the analogy with strong and fragile glasses [1], we can argue that the decrease of T 0 with increasing H signals that the apparent glassy state (at T < T 0 ) becomes stronger with H. This is consistent with the magnetic field increasing the effective disorder, which further enhances the suppression of vortex diffusivity. The extrapolation of T 0 (H) to zero at H * 9 T suggests the existence of a quantum critical point separating the dissipationless solid (for H < H * ) from the SC vortex glass phase that only exists at T = 0. Interestingly, a similar field-tuned transition between two superconducting ground states with different ordering of the vortex matter, a vortex solid at lower H and a T = 0 vortex glass at higher H, has been observed also in underdoped copper-oxide high-temperature superconductors [37,38], which are relatively clean quasi-2D systems. Further insight into the physical mechanisms responsible for the fragile-glass dynamics of the melting of a 2D VL could come from experiments on more disordered thin films and from theoretical exploration of the evolution of the dynamical heterogeneities in the presence of finite disorder.

METHODS
Samples. a-MoGe films with thickness t = 22 nm were grown on surface-oxidized Si substrate through pulsed laser deposition. The reported chemical stoichiometry of these thin films, as seen from the dispersive X-ray analysis, is Mo 71±1.5 Ge 29±1.5 ; their properties have been described in detail elsewhere [8,9,18]. The films were capped with a 2 nm thick Si layer to prevent surface oxidation, and patterned in Hall bar geometry using a shadow mask. Detailed measurements were performed on two samples with dimensions 0.36 mm (width)×2.8 mm (length), and 1.1 mm distance between the voltage contacts. The two samples exhibited an almost identical behavior. For sample 1, the voltage contact width is 0.05 mm and the zero-field T c = (7.70 ± 0.05) K; for sample 2, the voltage contact width is 0.025 mm and T c = (7.50 ± 0.05) K. T c is defined as the temperature at which the resistivity start to rise above the experimental noise floor (∼ 3.6 × 10 −4 mΩcm or ∼ 0.5 Ω in resistance). Gold leads (≈ 50 µm in diameter) were attached to the samples (on top of a Si layer) using the two-component EPO-TEK-E4110 epoxy. The resulting contact resistances were ∼ 200 Ω each at room temperature.
Measurements. Resistance was measured using the standard four-probe ac technique (∼13 Hz or 17 Hz) with either SR 7265 lock-in amplifiers or a LS 372 resistance bridge. Some of the measurements were performed using a dc reversal method with a Keithley 6221 current source and a Keithley 2182A nanovoltmeter. The excitation current densities were 0.1-10 A cm −2 , depending on the temperature, and low enough to avoid Joule heating. dV /dI measurements were carried out by applying a dc current bias I dc and a small ac current excitation (∼ 13 Hz) through the sample (I ac = 1 µA at T > 1 K and I ac = 10 nA at T < 1 K), while measuring the ac voltage across the sample for 150 s and recording the average value for each I dc .
Several cryostats were used to cover a wide range of temperatures and fields: a HelioxVL 3 He system (0.25 ≤ T ≤ 200 K) with H up to 9 T; a dilution refrigerator (0.02 ≤ T 1 K) and a 3 He system (0.3 ≤ T < 60 K) in superconducting magnets with H up to 18 T; and a variable-temperature insert (1.8 < T ≤ 200 K) in a Quantum Design PPMS with H up to 9 T. The fields, applied perpendicular to the film surface, were swept at constant temperatures. A low sweep rate of 0.1 T/min was used to avoid heating of the sample due to eddy currents. Measurements in the dilution refrigerator and the HelioxVL 3 He system were equipped with filters, consisting of a 1 kΩ resistor in series with a π-filter [5 dB (60 dB) EMI reduction at 10 MHz (1 GHz)] in each wire at the room temperature end of the cryostat to reduce high-frequency noise and heating by radiation. Furthermore, the dilution refrigerator and the 3 He system measurements in 18 T magnets were carried out in the Millikelvin Facility of the National High Magnetic Field Laboratory, which is an electromagnetically shielded room.
Simulations. We have performed Monte Carlo (MC) simulations on a spin system on a square grid with lattice spacing a = 1, linear size L = 56 and uniform magnetic field, whose intensity can be expressed in terms of the quantum-flux fraction f passing through a unitary plaquette f = Ba 2 /Φ 0 .
Here we considered f = 1/L, which results in N v = f L 2 = 56 vortices with a given vorticity. Although the system size simulated is smaller than in other numerical simulations of particle systems [10,[32][33][34][35][36], it is the state of the art for numerical simulations of vortex lattices within XY models. Each MC step consists of the local updating of all spins of the lattice by means of the Metropolis-Hastings algorithm. All observables have been computed at equilibrium, as achieved after a certain numbert of MC steps. For temperatures above T 0 , we identifȳ t with the entrance into the diffusive regime, while below T 0 we estimated thatt 10 8 MC steps provides a stable result. After discarding the firstt steps, we reset the Monte Carlo time and proceed with the measurements of both the statical and the dynamical observables.
The superfluid stiffness J s is defined as the linear response to an infinitesimal phase twist in a given direction, say µ, and it reads: It accounts for two contributions: the diamagnetic part J µ d , proportional to the energy density of the system, and the paramagnetic part J µ p , given by the connected current-current response function. As in Eq. (5), here and in what follows, . . . stands for both the thermal average, performed at equilibrium over all the MC steps, and for the average over ten independent numerical simulations.
The local orientational order parameter ψ 6j is obtained by means of a Delaunay triangulation of the VL. It is defined for the hexagonal symmetry, as where N j is the number of nearest neighbors of the j-th vortex, and θ jk is the angle that the bond connecting the two neighboring vortices j and k forms with respect to a fixed direction in the plane. The global six-fold orientational order parameter G 6 is then obtained by summing over all the N v vortices and by computing its average, as given by Eq. (5). The susceptibility of the orientational oder parameter χ 6 is defined as where Ψ 6 = 1 Nv Nv j=1 ψ 6j with ψ 6j defined in (9). The temperature dependence of χ 6 ( Fig. 11) exhibits a peak at T = 0.05 that we identify as the critical temperature T hex between the hexatic and the isotropic liquid phase. In both panels, the Arrhenius (red dashed lines) and VFT (blue dotted lines) fits are performed to the data shown within the two red arrows. The VFT form describes the data better, especially at the lowest T . SinceT0 decreases with increasing H, the difference between the two fits also decreases at higher fields. For the same data, fits to the diffusion coefficient, Eq. (3), yield T0 = (4.69 ± 0.02) K and T0 = (1.14 ± 0.06) K for 3 T and 8.2 T, respectively; these values are the same within error as the correspondingT0.
FIG. 6. The VFT fitting parameters T0, W , and ln Z. The data are shown for sample 1. The error bars represent standard errors obtained from the nonlinear fits of ( 2e h ) 2 kBρT to the VFT law in Eq. (3). A strong correlation between W (purple triangles) and ln Z (magenta squares) is observed (ln Z is multiplied by 2.5 on the y-axis for clarity). T0 (blue dots) is independent of the other two parameters. 7. Comparison of the parameters W (from the VFT fit) and U (from the Arrhenius fit). The data are shown for sample 1. W (purple triangles) and U (black diamonds) differ considerably for H < 9 T, but they become equal, within error, for H 9 T (see inset) where T0 = 0, thus confirming the consistency of the analysis. Dashed lines guide the eye. The error bars for W and U are standard errors from the nonlinear fit of ( 2e h ) 2 kBρT to the VFT form and linearized fit of ρ to the Arrhenius form, respectively. FIG. 9. Differential resistance dV /dI as a function of dc current I dc in different regimes. a, In the "Arrhenius" region ( Fig. 2) at T = 1.3 K and H = 9 T, a strong nonlinearity is observed, consistent with the motion of vortices in the presence of disorder; Tc(H = 9 T) = 1.1 ± 0.2 K. b, Similar nonlinear dV /dI vs I dc in the "Arrhenius" region for H = 12.7 T and T = 0.017 K. At a much higher field (H = 18.0 T> Hc2), in the normal state, dV /dI vs I dc is ohmic. The nonlinear behavior becomes unobservable also as T increases, as shown in c, for H = 9 T at T = 1.5 K, and d, for H = 9 T at T = 1.75 K, at lower currents than in a. e, The evolution of dV /dI vs I dc with increasing H, as shown, at T = 0.017 K. The data for H = 18.0 T and H = 12.7 T are the same as in d. When ρ drops below the noise floor, e.g. for H = 12.25 T, a finite I dc 2 µA, i.e. a critical current, is needed to depin a measurable number of vortices. At somewhat higher H = 12.5 T, ρ is nonzero at low current bias, consistent with the vortex creep due to thermal fluctuations. The depinning effect becomes observable at the critical current I dc ∼ 0.1 µA, which is lower than that for H = 12.25 T, as expected. f, dV /dI vs I dc is ohmic at T ≈ 5.5 K and H = 2.5 T, where ρ(T ) obeys the VFT law. The error bars in all panels correspond to ±1 SD obtained from averaging the ac voltage over 150 s at a fixed I dc . Dashed lines guide the eye. The T fluctuations for these measurements were less than 5 mK.

S1. MONTE CARLO SIMULATIONS: STATIC CHARACTERIZATION OF THE VORTEX LATTICE
To further characterize the three phases found, we have computed the structure factor of the vortex lattice defined as: where N v is the total number of vortices, ν(r i ) is the local vortex density, equal to one if a vortex occupies the site r i and zero otherwise, and k is the vortex lattice reciprocal vector. In Fig. 10, we show S(k) computed at three different temperatures: T = 0.02, 0.04, 0.07, corresponding respectively to the solid, hexatic-liquid, and isotropic-liquid phase. At high temperature, the structure factor presents the typical circular symmetry of an isotropic liquid. Decreasing the temperature, such symmetry breaks down and the Bragg-peaks structure appears showing six well-defined spots. Finally, in the solid state the Bragg peaks become well-defined even at large k. Let us highlight that the structure-factor anisotropy observed at low temperature is due to the commensurability of the VL with respect to the underlying square lattice of spins. The VL has indeed two ways to align with respect to the underlying square lattice: either to be perfectly commensurate with thex-axis or with theŷ-axis. In order to identify the isotropic to hexatic liquid transition, we have computed the orientational susceptibility χ 6 . The susceptibility of the orientational oder parameter χ 6 is defined as where Ψ 6 = 1

Nv
Nv j=1 ψ 6j with ψ 6j defined in Eq. (9) of the main text. The temperature dependence of χ 6 (Fig. 11) exhibits a peak at T = 0.05 that we identify as the critical temperature T hex between the hexatic and the isotropic liquid phase.

S2. MONTE CARLO SIMULATIONS: SIGNATURES OF HETEROGENEOUS DYNAMICS
The heterogeneous nature of the vortex dynamics in the hexatic phase has its fingerprints in different observables. Together with the mean-square displacement ∆r 2 (t) , reported in the main text, we also computed the self-part of the intermediate scattering function F s (|k * |, t), and the non-Gaussian parameter α 2 (t). The resulting trends in time for different temperatures are shown in Fig. 12, where the time variable t labels the discrete MC steps. To highlight the onset of the hexatic phase, the curves at the verge of the isotropic-liquid to the hexatic-liquid transition (T = 0.05) have been plotted in gray. The self-part of the intermediate scattering function is a dynamical autocorrelation function for the VL and it is defined as: where |k * | = 2π/a v (with a v the lattice spacing of the VL) is the reciprocal vector at which the structure factor shows its first peak (see Fig. 10), r j (t) is the position of the j-th vortex at the MC time t and t M is the largest MC time used in the simulations. The thermal average is thus the sum over all the possible t 0 for a given time t, while (. . . ) stands for the average over ten independent simulations. From Fig. 12a, one can see that, in the isotropic liquid phase, F s (|k max |, t) decays exponentially to zero with a unique relaxation time. On the other hand, as the hexatic phase is approached, it starts showing a plateau, which increases with decreasing T . This two-step relaxation decay is another typical signature of a caging mechanism [1] due in this case to the onset of the hexatic phase. The mean-square displacement, already shown in the  inset of Fig. 3b, is defined as: As discussed in the main text, at the onset of the hexatic liquid phase an intermediate subdiffusive regime appears (Fig. 12b). This signals an inhibition of the particle motion due to the cage formed by the neighboring particles, similarly to what happens in supercooled liquid and glassy systems. The intermediate plateau observed both in the mean-square displacement and in the intermediate scattering function signals the emergence of a heterogeneous dynamics within a certain time scale. To further investigate this behavior, we have computed the non-Gaussian parameter α 2 (t) = 1 2 ∆r 4 (t) which quantifies the heterogeneity of the dynamics in terms of strength and extent in time [2]. We find (Fig. 12c) that, at very short times, the VL displays a strongly heterogeneous dynamics for all the temperatures analyzed. This anomalous behavior is due to the underlying squared grid of spins that forces vortices to move only along four possible directions: ±x and ±ŷ. With increasing time, the direction of motion becomes less sensitive to the underlying grid and α 2 (t) decreases. Apart from the initial heterogeneity, at high temperatures the VL shows a homogenous dynamics. On the other hand, at the onset of the hexatic phase, the non-Gaussian parameter starts displaying a dome at longer time scales. This additional signature of heterogeneity is another indication of a caging mechanism triggered by the onset of the hexatic phase.
With decreasing temperature the height of the dome increases, signalling the increase of the heterogeneous dynamics strength. At the same time, the peak of the dome moves to longer times, reflecting the increase of the time scale at which vortices escape from their cage.

S3. BRIEF DESCRIPTION OF THE SUPPLEMENTARY VIDEO
Supplementary Video S1: Heterogeneous vortex motion in the hexatic phase. Images of the vortex motion at T = 0.045 for different MC times. In this regime, the dynamics is strongly heterogeneous in both space and time. The vortices remain indeed trapped for a long time in the cage formed by their neighbors before exiting in a collective burst along the symmetry axis of the vortex lattice.
In the video, the vortices are represented as colored disks, centered at a given position x i at time t i . Each frame shows vortex evolution over ten consecutive discrete time steps t i . To help visualize the time evolution, we assigned to each disk a radius r i that is increasing with increasing time t i . As a consequence, larger disks identify the vortex position at larger times. In addition, we add a solid black line connecting the centers of the disks to help visualize the vortex motion as a function of time. The grey lines in the background show the trails left by the vortices during the whole simulation. The horizontal line at the bottom of the image indicates the time flow. Note that the time steps t i are not equally spaced with respect to the MC time. Indeed, for the sake of the memory allocation, we have stored the data each t MC = int(A q ) + kA Nq , with q ∈ [0, N q ], k ∈ [0, N k ] and A = 1.3.