Dynamic scaling in the 2D Ising spin glass with Gaussian couplings

TL;DR

This study uses simulated annealing and Kibble-Zurek scaling to analyze the dynamic critical behavior of the 2D Gaussian-coupled Ising spin glass, revealing a large dynamic exponent and differences from bimodal coupling systems.

Contribution

It provides the first detailed dynamic scaling analysis of the 2D Gaussian Ising spin glass, highlighting the impact of coupling distribution on relaxation dynamics.

Findings

Dynamic critical exponent z=13.6 ± 0.4 for both order parameter and energy.

Different relaxation behaviors compared to bimodal couplings due to ground state degeneracy.

Kibble-Zurek scaling applies, but with altered scaling functions due to quasi-degenerate states.

Abstract

We carry out simulated annealing and employ a generalized Kibble-Zurek scaling hypothesis to study the 2D Ising spin glass with normal-distributed couplings. The system has an equilibrium glass transition at temperature $T=0$. From a scaling analysis when $T\rightarrow 0$ at different annealing velocities, we extract the dynamic critical exponent $z$, i.e., the exponent relating the relaxation time $\tau$ to the system length $L$; $\tau\sim L^z$. We find $z=13.6 \pm 0.4$ for both the Edwards-Anderson spin-glass order parameter and the excess energy. This is different from a previous study of the system with bimodal couplings [S. J. Rubin, N. Xu, and A. W. Sandvik, Phys. Rev. E {\bf 95}, 052133 (2017)] where the dynamics is faster and the above two quantities relax with different exponents (and that of the energy is larger). We here argue that the different behaviors arise as a consequence of the different low-energy landscapes---for normal-distributed couplings the ground state is unique (up to a spin reflection) while the system with bimodal couplings is massively degenerate. Our results reinforce the conclusion of anomalous entropy-driven relaxation behavior in the bimodal Ising glass. In the case of a continuous coupling distribution, our results presented here indicate that, although Kibble-Zurek scaling holds, the perturbative behavior normally applying in the slow limit breaks down, likely due to quasi-degenerate states, and the scaling function takes a different form.

Equations

$$v_{\mathrm{KZ}}\propto L^{-z-1/\nu },$$

$$H=\sum _{⟨ij⟩}{J}_{ij}{\sigma }_i{\sigma }_j,{\sigma }_i=\pm 1,$$

$$q=\frac{1}{N}\sum _{i=1}^N{\sigma }_i^{(1)}{\sigma }_i^{(2)},$$

$$E(L)-E_{\infty }=a{L}^{-(d+\frac{1}{\nu })}.$$

$$T(t)=T_{\mathrm{ini}}(1-t/t_{\max }{)}^r.$$

$$A(v,L)=A_{\mathrm{eq}}(L)f(v/v_{\mathrm{KZ}}),$$

$$⟨q^2(v,L)⟩\propto \{\begin{array}{l}{f}_0(v{L}^{z+1/\nu }),v\lesssim v_{\mathrm{KZ}},\\ \\ (v{L}^{z+1/\nu }{)}^{-x},v_{\mathrm{KZ}}\lesssim v\lesssim 1,\\ \\ L^{-2}{f}_1(1/v),v\gtrsim 1.\end{array}$$

$$x=\frac{2}{z+1/\nu },$$

$$\Delta E(v,L)=E(v,L)-E_{\infty },$$

$$⟨\Delta E(v,L)⟩\propto \{\begin{array}{l}{L}^{-(2+1/\nu )}{g}_0(v{L}^{z+1/\nu }),v\lesssim v_{\mathrm{KZ}},\\ \\ L^{-(2+1/\nu )}(v{L}^{z+1/\nu }{)}^{-x^{\prime }},v_{\mathrm{KZ}}\lesssim v\lesssim 1,\\ \\ g_1(1/v),v\gtrsim 1,\end{array}$$

$$x^{\prime }=\frac{2+1/\nu }{z+1/\nu }.$$

$$⟨q^2⟩=1-bL/{\xi }_v,(L/{\xi }_v\to 0),$$

$${\xi }_v\propto v^{-1/(z+1/\nu )},$$

$$E(v,L)$$

Full text

I Introduction

Spin glasses are benchmark models for studying complex physical systems and optimization problems. Due to the disorder and frustration (random mixed-sign couplings), the energy landscapes of these systems are very rough, with many local minimums, and it is very challenging to find the true global minimum (ground state) through Monte Carlo (MC) simulations edwards75 ; binder86 ; fisher91 ; parisi83 . Among the common spin glass systems, the two-dimensional (2D) Ising spin glass (2DISG) is special in that the paramagnetic–glass phase transition occurs exactly at temperature $T=0$. The system has long-range spin-glass order (defined with the Edwards-Anderson, EA, order parameter) at $T=0$, and the correlation length diverges as a power law, $\xi\sim T^{-\nu}$ when $T\rightarrow 0$. Many works have been devoted to the nature of the critical behavior and to obtaining the critical exponents of 2DISG system with both normal-distributed (Gaussian) and bimodal couplings kawashima92 ; campbell04 ; houdayer04 ; rieger96 . However, due to the considerable challenges with MC simulations, especially for large systems at low temperature, there are still significant issues under debate. For example, whether or not the 2DISG with bimodal $J=\pm 1$ and Gaussian couplings belong to the same universality class in their equilibrium criticality is still in question jorg06 ; toldin11 ; thomas11 ; fernandez16 ; campbell16 ; campbell17 . Undisputed is the fact that the ground-state properties of the two models are different. The system with Gaussian couplings has a unique (non-degenerate) ground state, up to a trivial spin reflection, while the model with bimodal couplings has infinite degeneracy in the thermodynamic limit.

Given the difficulties in studying the critical behavior through equilibrium simulations, the recently developed non-equilibrium approach based on generalized Kibble-Zurek (KZ) scaling kibble76 ; zurek85 ; zhong05 ; chandran12 ; liu14 ; liu15a ; polkovnikov05 ; zurek05 provides a powerful alternative method for studies of spin-glass models. KZ scaling of simulated annealing (SA) results has been successfully applied to 3D and 2D spin glass systems in order to extract the dynamic exponent $z$ and other critical exponents liu15 ; rubin17 . The key aspect of the KZ mechanism used in this context is the prediction that slow (close to equilibrium) and fast (far from equilibrium) SA processes are separated by an annealing velocity $v_{\rm KZ}$ that scales with the system size (length) $L$ as

[TABLE]

where $\nu$ is the standard equilibrium correlation-length exponent. Here we apply this approach to the 2DISG with Gaussian couplings, following the recent work on bimodal couplings rubin17 .

In Ref. rubin17, , a surprising behavior with dual time scales governing the relaxation when $T\to 0$ was discovered. Contrary to the general expectation that the order parameter is the slowest-relaxing physical observable, and that most other quantities are asymptotically governed by that same time scale, a larger dynamic exponent, $z_{E}\approx 10.3$, was found for the excess energy than $z_{q}\approx 8.3$ for the EA order parameter. The physical mechanism proposed to underly the two time scales relies on the backbone (largest common cluster) and droplet (zero-energy flippable cluster) structure of the massively degenerate ground states of the $J=\pm 1$ model roma10 ; thomas11 , which leads to a concentration in the configuration space of low-energy states that entropically attracts the SA process. The proximity of true ground states and low-energy excitations to each other within this region was proposed to lead to an insensitivity of the replica-overlap definition of the order parameter to low-energy excitations, so that the final relaxation of the energy leads to only sub-leading corrections to the already equilibrated mean order parameter. This $T\to 0$ relaxation process is of particular relevance in related optimization problems, where currently there is much interest in comparing SA and quantum annealing protocols and the connectivity (especially the dimensionality) of the the spins (qubits in quantum annealing) may play a very important role katzgraber15 .

It should be noted that the relaxation dynamics in an SA process for $T\to 0$ can be very different from the dynamics associated with ergodic sampling at fixed $T>0$. The latter should be associated with a divergent dynamic exponent when $T\to 0$ in 2D Ising spin glasses katzgraber05 , which also is consistent with the non-ergodicity of local spin moves at $T=0$. In SA, the temperature is constantly changing and naive arguments based on activated dynamic scaling to overcome energy barriers do not necessarily apply in all cases, since the details of the energy landscape matter. In Ref. rubin17, it was argued that the droplet structure of the 2D Ising spin glass corresponds to a funnel-like feature of the energy landscape where high energy barriers can be overcome at high temperatures and the barriers remaining as the ground state is approached when $T\to 0$ become typically smaller, such that a power-law scaling of the annealing time required to reach the ground state is obtained. This situation is also of great interest in the context of optimization and computational complexity, as a case where the typical exponential scaling to find an optimal solution can be avoided katzgraber15 .

In the case of Gaussian-distributed couplings, which we study in this paper, the backbone structure can be defined only as an approximation with low-energy states instead of true ground states roma13 ; roma16 . Strictly speaking, there is no definable backbone and zero-energy clusters in that model due to the lack of ground-state degeneracy. Because of this qualitative difference of the ground-state landscape, one can expect different dynamical properties for the Gaussian model (or any other continuous coupling distribution). The aim of the work presented here is to apply exactly the same scaling approach as was done with the bimodal 2DISG in Ref. rubin17, and test whether a clearly different asymptotic relaxation mechanism can be detected. We will show that, indeed, in this case the excess energy and the EA order parameter relax with the same dynamic exponent, and the value of the exponent, $z=13.6(2)$ (where here and later the number in parentheses indicates the statistical error of the preceding digit), is significantly larger than both exponents found in the bimodal case.

The organization of the rest of the paper is as follows: In Sec. II we discuss the known equilibrium properties and expected finite-size behaviors near the $T=0$ critical point of the 2DISG model with Gaussian couplings. These properties are important when extending finite-size scaling to non-equilibrium setups where the annealing velocity enters as another variable. We describe the SA simulation procedures, where we have applied GPU (graphics processing-unit) computing for very efficient MC sampling with the Metropolis algorithm, and summarize the KZ scaling procedures we have applied to quantify the relaxation behavior as a function of system size and annealing velocity. In Sec. III we present results of the scaling analysis for the excess energy and the EA order parameter. Lastly, in Sec. IV we further discuss our findings and contrast them with the conclusions previously drawn for the bimodal case.

II Model and methods

The Hamiltonian of the 2DISG is

[TABLE]

where, in the case considered here, $\langle ij\rangle$ stands for nearest-neighbour spins on a 2D square lattice with $L^{2}$ sites and periodic boundary conditions. The couplings $J_{ij}$ are drawn from some distribution, here Gaussian with mean [math] and standard deviation $1$.

II.1 Equilibrium finite-size scaling

The primary quantity capturing the spin-glass phase transition is the EA order parameter,

[TABLE]

where (1) and (2) stand for two independently generated configurations (two different MC simulations), also referred to as ‘replicas’, of systems with the same coupling realization $\{J_{ij}\}$. In this paper, we focus on the mean squared EA order parameter, $\langle q^{2}\rangle$, as well as the internal energy density $E=\langle H\rangle/N$ in the limit $T\to 0$ reached in SA simulations with Metropolis dynamics. For simplicity of notation, we use $\langle...\rangle$ to denote the combined MC expectation value and the average over disorder samples.

It is known that the 2DISG with Gaussian couplings has a phase transition exactly at $T=0$, and its critical behavior has been studied extensively kawashima92 ; campbell04 ; houdayer04 . Unlike the 2DISG with $J=\pm 1$ couplings, where there are many degenerate ground states, there is only a unique ground state (and the state with all spins reversed). Thus, as $T\rightarrow 0$ all the independent replicas will eventually fall into the same ground state configuration in the limit of a very slow SA process, and the EA order parameter $\langle q^{2}\rangle$ must approach $1$ without any finite-size corrections in the $T=0$ value. However, according to the study in campbell04 , the equilibrium ground-state energy density has a finite-size correction of the form

[TABLE]

Here the energy per spin for infinite $d=2$ dimensional system is $E_{\infty}=-1.314788(4)$ thomas07 and the most precise value available for the critical exponent $\nu$ of the correlation length $\xi$ (where in the case $T_{c}=0$ we have $\xi\sim T^{-\nu}$) is $\nu=3.56(2)$ campbell04 . The prefactor $a$ of the scaling in $L$ was claimed to be exactly $a=1$. In the following analysis of SA data, we will make use of the form (4) with the previously determined values of $E_{\infty}$ and $\nu$ (while the value of $a$ is less important).

II.2 Simulated annealing

Most of the simulations were run on Nvidia CUDA enabled GPUs, with single-spin Metropolis updates and multi-spin coding where the Ising spins $\sigma_{i}=\pm 1$ of the model [Eq.(2)] are coded as bits of 32-bit integers. Thus, with the same set of random couplings, one simulation propagates $32$ replicas from different initial conditions and the order parameter $q^{2}$ is computed at the end of the run (at $T=0$) from the overlap, of the form Eq. (3), among these replicas. One sweep of MC updates involves $N=L^{2}$ Metropolis spin-flip attempts, carried out successively in two groups corresponding to the standard checker-board decomposition of the lattice, so that all spins in a given sublattice can be updated in parallel independently of each other. For each of the 32 replicas, different random numbers are generated in order for the updating processes to be fully independent. Our program achieves around $2\times 10^{9}$ attempted spin flips per second on a single GPU. Further discussion on how to implement these updates on GPUs can be found in Refs. preis09, ; block10, ; lulli15, ; hsieh13, ; hartmann15, ; fang14, ; jesi14, ; parisi10, ; weigel12a, ; weigel12b, .

In an SA procedure, after initial equilibration of the system at a high temperature $T_{\rm ini}$, the temperature $T(t)$ is lowered as a function of the simulation time $t$ according to some protocol. In the context of the KZ mechanism one normally considers the approach to a phase transition using a linear protocol, or, if the transition point is known one can approach it with a nonlinear power-law protocol (or in principle some other protocol). Note that we are here not interested in finding an optimal SA protocol (i.e., the one that would bring us to the ground state in the shortest time), but aim to test the power-law KZ scaling hypothesis for the 2D Ising glass with transition temperature exactly at $T=0$, and use it to extract dynamic information.

The general power-law protocol takes us to $T=0$ as a function of the total simulation time $t_{\rm max}$ according to

[TABLE]

Here $r=1$ corresponds to the standard SA protocol where the temperature is decreased linearly. In order to disentangle the exponents $z$ and $\nu$ involved in KZ scaling, e.g., in Eq. (1), it is also useful to study other values of $r$, as exemplified in the previous study of the 2DISG with bimodal couplings rubin17 . There a consistency check was provided by the fact that the entropy exponent $\Theta_{S}$, which plays the role of $1/\nu$ in that case thomas11 , was determined independently and agreed with previous calculations. In the work reported here, we only consider $r=1$ and use the known value of $\nu$ to extract $z$, because the calculations with Gaussian couplings are very expensive (even with the use of GPUs), mainly due to the fact that longer times are needed to reach close to the unique ground state. Another reason for only considering $r=1$ is that the value of $1/\nu$ is small, around $0.28$, and it is then hard to determine it independently from simulations of two or more $r$ values within the error levels we can reach for the KZ exponent, $z+1/\nu$ in Eq. (1) for $r=1$ [and $z+1/r\nu$ for other $r>0$ polkovnikov05 , which we do not consider here].

The annealing velocity is defined as $v=T_{\rm ini}/t_{\rm max}$. At the last step of the annealing process, when $T$ has reached [math], we take measurements of the EA order parameter $q^{2}$, the energy (per site) $E$ of the system, as well as the minimum energy per spin, $E_{\rm min}$, among any of the $32$ replicas. The SA process is repeated many times for different realizations of the random couplings.

To test the SA program, for $L\leq 6$ we used sufficiently long simulations for several disorder realizations to relax the systems all the way to the ground state. We checked these results against exact ground states, which can be obtained using an exhaustive search in the state space or using a matching algorithm such as those described in Refs. thomas07, ; hartmann11, . Based on the tests we know that the ground states are indeed reached in the simulations for sufficiently slow $v$, as expected. For the mean values taken over a large number of samples that we report below, because of the long simulation times we were not able to use low enough $v$ for even the small-$L$ systems to reach their ground state in all cases, and for larger $L$ none of the systems reached the ground state. As we will see below, the mean values still do reach sufficiently close to their ground-state values to test the asymptotic KZ relaxation behavior.

II.3 Dynamic scaling

In a generalized KZ scaling ansatz, for a system reaching the critical point through the annealing protocol expressed in Eq. (5), a physical quantity $A$ evaluated at the critical point can be written in the following finite-size scaling form zhong05 ; chandran12 ; liu14 :

[TABLE]

where the “critical” KZ velocity for the linear SA protocol, $r=1$ in Eq. (5), is given by Eq. (1) up to an undetermined and essentially arbitrary factor. This velocity demarks the borderline between fast and slow annealing processes. The function $A_{\rm eq}(L)$ in Eq. (6) stands for the equilibrium finite-size dependent quantity $A$ at the critical point, which normally is a power of $L$ to leading order but also can include scaling corrections. The dynamic exponent relates the scaling of the relaxation time to the correlation length through $\tau\sim\xi^{z}$, which at the critical point for finite-size systems turns into $\tau\sim L^{z}$ by the standard substitution $\xi\to L$ in critical finite-size scaling. Recall the discussions in Sec. II.1, that when the system is in equilibrium at $T=0$ the order parameter $\langle q^{2}\rangle=1$ without any finite-size effect, while the excess energy density has a finite-size correction of the form Eq. (4). These behaviors will be reflected in the corresponding $A_{\rm eq}(L)$ in Eq. (6).

According to the general non-equilibrium scaling form that describes the dynamics in its full regime of velocities and sufficiently large system sizes, the order parameter $\langle q^{2}\rangle$ can be written in the following way liu14 :

[TABLE]

Here the first line describes the slow velocity regime, where the function $f_{0}$ normally would be a regular (Taylor-expandable) function of the KZ-scaled velocity $vL^{z+1/\nu}$, although below we will argue that, in the case considered here, a corresponding function of a power of the KZ variable has to be used for this to be true. As discussed above, there should not be any dependence on the system size asymptotically for $v\to 0$ since $\langle q^{2}\rangle\to 1$ because of the unique ground state. In principle there could be $L$-dependent corrections for $v>0$ but the form of these are not presently known. The third line describes the fast velocity regime, in which the system size is larger than the correlation length $\xi_{v}$ at the end of the annealing process and, thus, there is no dependence on $L$, other than the trivial factor $L^{-2}$ that follows from Eq. (3) when the spin-glass correlation length is finite. The function $f_{1}$ should be Taylor-expandable in $1/v$. The second line in Eq. (7) describes the intermediate power-law regime that connects the two other regimes. It follows from the scaling hypothesis, Eqs. (6) and (1), where the behavior when $L\to\infty$ at fixed $v$ must reduce to (connect smoothly to) the form on the third line, again because $\xi_{v}\ll L$ in this limit. The only way to make this possible (i.e., to ensure the same $L$-dependence in the two forms) is with the power-law form $f(v/v_{\rm KZ})\to$ $(vL^{z+{1}/{\nu}})^{-x}$, where the exponent $x$ must be given by

[TABLE]

so that the power-law form can also be written as $\langle q^{2}\rangle\propto L^{-2}v^{-x}$. Then the connection between lines 2 and 3 in Eq. (7) corresponds to the function $f_{1}(1/v)$ crossing over into the form $v^{-x}$ and the connection between lines 1 and 2 corresponds to $f_{0}(vL^{z+1/\nu})$ taking the form $(vL^{z+{1}/{\nu}})^{-x}$ for large $vL^{z+{1}/{\nu}}$. In other words, the KZ scaling form in Eq. (6), with the KZ velocity given by Eq. (1), covers the first and second lines of Eq. (7), while the third line represents the break-down of this form for higher velocities.

We also consider the excess energy density, which we here define relative to the known infinite-size equilibrium $T=0$ value $E_{\infty}$,

[TABLE]

i.e., it contains contributions from both finite size and non-zero velocity. In analogy with the above discussion of the EA order parameter, and considering the equilibrium finite-size scaling given in Eq. (4), the behaviors in the three different velocity regimes should be given by

[TABLE]

where, unlike Eq. (10), there is no $L$ dependence on the third line because the excess energy is defined per spin and takes a constant value when $v\to\infty$ (i.e., in the initial state). In this case, for the power-law regime to be valid, i.e., for there to be no size dependance on the second line ($\Delta E\sim v^{-x^{\prime}}$), the exponent $x^{\prime}$ is given by

[TABLE]

In the next section, we will present our results of the application of the above scaling forms.

III Results

All simulations reported here started from $T_{\rm ini}=8$, where the system can be easily equilibrated. Starting from a random configuration for each disorder sample, we used 10 MC sweeps at this initial temperature. From there, we used the linear SA process, i.e., $r=1$ in Eq. (5), and measurements were taken at the last step of the annealing process where $T=0$. We used system sizes from $L=4$ to $L=64$. To span a wide range of velocities, we take the total time for the simulations as $t_{\rm max}=2^{n}$, where $n=2,3,...,30$ for small system sizes, while for large system sizes we only used $n$ up to $28$ to stay within reasonable computing times. To obtain good statistical averages, we simulated at least $5\times 10^{3}$ coupling realizations in most cases and $10^{3}$ realizations for the lowest velocities and largest system sizes.

III.1 Mean excess energy density

Figure 1(a) shows the velocity scaling of the average of the excess energy density, Eq. (9), with $E_{\infty}=-1.31479$ from Refs. campbell04, ; thomas07, . The overall expected size dependence in equilibrium from Eq. (4) has been divided out, and the velocity has been rescaled according to the expected KZ form in Eqs. (6) and (1). Here we use data points from system sizes $L=8$ to $L=64$ in the data-collapse procedure, for each $L$ excluding velocities too high to give results on the common scaling function. We vary the scaling exponent $z+1/\nu$ to achieve optimal collapse relative to a fitted polynomial, repeating the procedure many times with Gaussian noise added to the data points in order to compute the statistical error. We obtain $z+1/\nu=13.9(4)$. Since $\nu\approx 3.56$ campbell04 , the dynamic exponent governing the excess energy is $z=13.6(4)$.

In Fig. 1(a) it is clear that data points for larger $v$ systematically peel off from the collapsed function and the region of data collapse in the rescaled variable is pushed further out to the right as $L$ increases. Figure 1(b) shows the same data graphed according to the third line of Eq. (10). The data now collapse well for high velocities, and instead the data for slower velocities peel off systematically from the common function as equilibrium is approached for each system size (i.e., the correlation length $\xi_{v}$ becomes of the order of the system size). In both Fig. 1(a) and Fig. 1(b), the straight lines indicate power-law behavior as described in the second line of Eq. (10) with the expected slopes, $x^{\prime}$ and $-x^{\prime}$, respectively, given by Eq. (11).

In order to test alternatives to the KZ scenario, we have also analyzed the data in other ways. One might naively expect that the SA relaxation of the system should involve an exponentially long time scale when $T\to 0$, given energy barriers that have to be overcome when the system becomes trapped in local energy minimums. The resulting activated scaling is reflected in a divergent equilibrium dynamic exponent $z_{\rm eq}(T)$ when $T\to 0$ katzgraber05 . In Fig. 1(c) we test for activated scaling of SA in the thermodynamic limit by graphing the same data as in Fig. 1(b) versus $\ln(1/v)$ instead of $1/v$, still using logarithmic scales on both axes. On this plot a linear dependence would imply $\Delta E\sim\ln^{-a}(1/v)$ with some positive exponent $a$, instead of the behavior $\Delta E\sim v^{x\prime}$ that we argued for above. We do not see any clear-cut linear behavior on the log-log plot, with more curvature in the system-size converged data for the lowest velocities than in the KZ-scaled data in Fig. 1(b). While one could perhaps argue that the data approach a straight line also here, we point out that the KZ form $\Delta E\sim v^{x\prime}$ with a small exponent $x\prime\approx 0.17$ will inevitably look similar to the form $\Delta E\sim\ln^{-a}(1/v)$ in a limited window of the argument $\ln(1/v)$, because a small power looks very similar to a logarithm. Thus, if in the window in question we have $\ln(1/v)\sim v^{-b}$ for some small value of the exponent $b$, then the KZ form will look like $\Delta E\sim\ln^{-x\prime/b}(1/v)$, so that the exponent $a$ above is roughly $x\prime/b$.

Note again that the KZ scaling demonstrated in panels (a) and (b) of Fig. 1 is not merely relying on the power-law scaling in the limit $L\to 0$ in a rather small window of velocities that we have achieved, but is mainly manifested in the generalized finite-size scaling form that applies also when equilibrium is reached for the smaller system sizes in panel (a). Importantly, there is full consistency of the asymptotic slope in panel (b) with the exponent $x^{\prime}$ defined in Eq. (11) with the value of $z$ that also describes the data collapse to the left of the power-law regime in Fig. 1(a), i.e., the KZ scaling hypothesis also describes the deviations from the infinite-size collapsed curve for the smaller system sizes ($L=8$ and $12$) in Fig. 1(b). In combination with the previous results for the bimodal coupling distribution in Ref. rubin17, , where the dynamic exponent is smaller and the KZ behavior can be seen even more clearly, we take these results as strong evidence of KZ scaling also with the normal-distributed couplings. In the following sections we will present further extensive quantitative support for this scenario.

III.2 Minimum excess energy density

In Figs. 2(a,b) we present the velocity scaling of the minimum energy, $\Delta E_{\rm min}$, defined for each disorder sample as the lowest energy reached at $T=0$ among any of the $32$ replicas run in parallel. We fix the exponent $z=13.6$ to be the same as that for the average energy shown in Fig. 1. We see that the scaling also works very well here. If we instead treat the exponent as a variable and optimize its value for the best data collapse, we obtain $z=13.5(5)$ in excellent agreement (within the error bars) with the one previously obtained. Thus, as expected, the two energies scale in the same way and the agreement also serves as a consistency check on the procedures. Note that, although the dynamic exponent is the same, the scaling functions are clearly different. In Fig. 2(b) we plot out the two scaling functions in the same graph by only showing the data points that fall clearly on the collapsed curve. Given how the quantities are measured, at a given velocity, the minimum energy reached is always lower than (or in some cases equal to) the average energy after the final MC step. Based on a rough estimation from the two curves, $\langle\Delta E_{\rm min}\rangle$ relaxes about $10^{4}$ times faster to the asymptotic minimum value than $\langle\Delta E\rangle$. However, for larger values of the scaled velocity, and for sufficiently large system sizes, we expect the two energies to converge to the same power-law behavior with the exponent given by Eq. (11), and we see indications of this convergence as well in Fig. 2(b). We can also see that our results for $\Delta E_{\rm min}$ are consistent with the prefactor $a=1$ in the equilibrium size dependence, Eq. (4), as the scaled quantity is close to $1$ in the low-velocity limit (though $a$ may also be marginally above $1$).

III.3 Order parameter

We next turn to the EA order parameter. Figures 3(a,b,c) show different aspects of the scaling of $\langle q^{2}\rangle$ with the velocity and the system size. In Fig. 3(a), $\langle q^{2}\rangle$ is graphed against the KZ-scaled velocity, using the same value of the dynamic exponent as was extracted above using the excess energy. Here we cannot reach as close to the equilibrium behavior as for the energy (especially the minimum energy), but the approach of $\langle q^{2}\rangle$ to $1$ is still obvious and the data for the smaller system sizes collapse very well in this regime, as shown more clearly in the inset of Fig. 3(a). The expected pure power-law behavior for large arguments $vL^{z+1/\nu}$ is not yet reached with the system sizes accessible here—the corrections to the power law as the equilibrium behavior is approached appear to be much larger than in the energy. The behavior is nevertheless consistent with an approach to the predicted asymptotic power-law scaling (indicated by the line in the figure). We also carried out the data collapse procedure with $z$ as a free parameter, using system sizes $L=8-24$ for which sufficient overlaps in the scaling variable exist so that the data-collapse procedure is well-defined. Figure 3(b) shows a clear minimim in the $\chi^{2}$-value of the fit versus the scaling exponent, in very good agreement with the best exponent obtained for the energy scaling in Fig. 1. A full error analysis gives $z=13.6(2)$, which is consistent with but statistically better than $z=13.6(4)$ from the energy. Thus, in contrast to the bimodal 2DISG, where a difference in dynamic exponents for the two quantities was found to be $z_{E}-z_{q}\approx 2$ ($z_{E}\approx 10.3$ and $z_{q}\approx 8.3$) rubin17 , in this case a single exponent governs the relaxation dynamics (as we had fully expected for this case where the ground state is unique).

In Fig. 3(c) we analyze the high-velocity limit of the order parameter, which eventually should cross over into the power-law regime. Recall that the collapse of data graphed versus the velocity (here the inverse velocity) at high velocities is trivial, merely reflecting the correlation length at the end of the SA process being much less than the system size (in the limit of $v\to\infty$ simply being the correlation length of the starting high-temperature equilibrium state), so that there is no size dependence. The initial state determines the details of the corresponding function $f_{1}(1/v)$ on the third line of Eq. (7) at high velocities, before the cross-over into the universal form written explicitly on the second line. Here again, we see a very slow approach to the pure power law, similar to the cross-over from the low-velocity side, and we can only say that the behavior is consistent with the expected behavior with $z\approx 13.6$.

To investigate the approach to equilibrium in more detail, in Fig. 4(a) we analyze the deviation $1-\langle q^{2}\rangle$ of the EA order parameter from the asymptotic size-independent equilibrium value $1$. Here again we see good data collapse setting in from the left side of the graph and extending further to the right with increasing system size. In the region where $1-\langle q^{2}\rangle$ is small, the behavior follows a power law with a small, non-integer exponent. Here one would normally expect an integer exponent, corresponding to an analytic function $f_{0}(v/v_{\rm KZ})=f_{0}(vL^{z+1/\nu})$ on the first line of Eq. (7). This has been observed in KZ scaling studies of non-random isolated quantum systems under Hamiltonian dynamics degrandi11 , for which the leading power laws for different quantities were also derived using adiabatic perturbation theory. Here the value of the exponent $a\approx 0.073$ in the power law $(L^{z+1/\nu})^{a}$ is very close to half of the value of the exponent $x$ in Eq. (8). Assuming that $a=x/2=(z+1/\nu)^{-1}$, we see that the asymptotic form is

[TABLE]

where $\xi_{v}$ is the KZ correlation length corresponding to finite velocity in the thermodynamic limit polkovnikov05 ; zurek05 ; chandran12 ;

[TABLE]

which can also be simply obtained from Eq. (1) by replacing $L$ by $\xi_{v}$. Thus, we conclude that, unlike other cases studied so far degrandi11 ; liu14 , here $f_{0}(vL^{z+1/\nu})$ is not Taylor-expandable but a corresponding function $\tilde{f}_{0}(L/\xi_{v})$ is. We do not have an explanation for this apparently different analytic form of the scaling function in this case, but empirically the evidence is compelling, as seen more directly in Fig. 4(b) where we plot the data against $L/\xi_{v}$ and compare with a power-law with exponent exactly $1$, i.e., testing the asymptotic form Eq. (12).

One might perhaps question the claim that the observed power-law behavior in Fig. 4 should reflect the true asymptotic form, given that the scaling variable $vL^{z+1/\nu}$ is still very large in this region, roughly in the range $10^{4}-10^{8}$ in the power-law region. However, the alternative scaling variable $L/\xi_{v}$ is much smaller, of the order $1$. Since a scaling variable is always determined only up to some essentially arbitrary factor, a more relevant measure of closeness to the asymptotic behavior should be the value of the quantity studied. Considering that $\langle q^{2}\rangle$ is as large as $0.8$, or, in other words, in two typical replicas $\approx 90\%$ of the spins are the same, and approximately the same fraction of the spins should then be in their ground-state configurations. We would then expect that the remaining relaxation of a dilute concentration of spins should already be governed by the asymptotic form, although we cannot completely exclude a cross-over into a different form still closer to equilibrium. As we will see below, we can push a bit further into the low-velocity regime by considering smaller system sizes.

In the above analysis of the EA order parameter, the smallest system size used in Figs. 3 and 4 was $L=8$. For smaller sizes we see behaviors that can be explained only with substantial scaling corrections included. Figure 5 focuses on the scaling of $1-\langle q^{2}\rangle$ for small system sizes, from $L=4$ to $L=16$. In Fig. 5(a), even though the $L\geq 8$ data collapse well in a region of slow velocities with standard KZ scaling and the same value of $z$ used above, the data for $L=4$ and $L=6$ clearly deviate substantially from a common scaling function. Staying within the subset of possible scaling corrections with no velocity dependence, we add a correction to the KZ argument $vL^{z+1/\nu}$ by multiplying it with $1-aL^{-b}$, with $a$ and $b$ optimized for the best data collapse (keeping $z$ at the previous value). With $a\approx 1.7$ and $b\approx 0.4$, the data collapse is very good on the left side, where also the power-law behavior found previously is substantially extended, with no detectable change in the exponent. This gives added support to the power-law form corresponding to Eq. (12) indeed being the asymptotic behavior.

We have also tried to analyze the asymptotic approach of the energy density to its equilibrium value. Here we can in principle use the KZ ansatz following from the known equilibrium finite-size scaling form Eq. (4) written in the following way:

[TABLE]

where $f(x)\to 1$ when $x=vL^{z+1/\nu}\to 0$ and $g(x)\to 0$ in this limit. Using the form of the equilibrium value, $E(0,L)=E_{\infty}+aL^{-(2+1/\nu)}$, with the parameters determined previously campbell04 , as mentioned below Eq. (4), we can analyze $(E(v,L)-E(0,L))L^{2+1/\nu}$. Within the standard scenario it should be a Taylor-expandable function $g(x)$ without constant term for small values of $x$. Unfortunately, here our results from Fig. 1 (from which we just need to subtract $1$ if the factor $a$ above really is exactly $1$, which is certainly consistent with our data in Fig. 2) are not good enough (the statistical errors are too large) to extract any meaningful behavior in the low-velocity limit. We can therefore at present not determine whether an integer power in $x$ obtains, or whether the leading behavior is instead an integer power of $L/\xi_{v}$ as in the case of $1-\langle q^{2}\rangle$.

IV Discussion

We have studied relaxation dynamics in the 2DISG model with Gaussian-distributed couplings by carrying out SA simulations in the $T\rightarrow 0$ limit, where the system in equilibrium goes through a phase transition into the glass state. Through performing scaling analysis according to the KZ hypothesis, we were able to extract the dynamical exponents associated with the excess energy $\langle\Delta E\rangle$ and the EA order parameter $\langle q^{2}\rangle$.

For the excess energy density, defined with respect to a previously determined value in the thermodynamic limit campbell04 , a data-collapse analysis yields $z=13.6(4)$, and the same kind of scaling procedure applied to the order parameter gives $z=13.6(2)$. Thus, there is a unique time scale governing the relaxation of both the order parameter and the excess energy. This in itself is not unexpected (as long as one accepts that the KZ mechanism applies), but it is interesting in light of the recent discovery of two substantially different dynamic exponents in the 2DISG with bimodal couplings rubin17 . The heuristic explanation provided for that behavior relied on the massive degeneracy of the ground state, which is lacking in the case of couplings drawn from a continuous distribution. The ground state degeneracy has consequences for the relaxation of the mean order parameter as defined using replica overlaps. Considering that we here used the exact same kind of scaling procedures, our results for the Gaussian couplings also lend further support to the anomalous behavior in the bimodal case and its explanation in terms of ground-state degeneracy.

The dynamic exponent we find here for the system with Gaussian couplings is significantly larger than the two different exponents for the bimodal couplings, where the larger of the two dynamic exponent, i.e., the one governing the energy relaxation, is $z^{\prime}\approx 10.4$. While we do not have a rigorous explanation for this difference, it should be related to the fact that the ground state in the case of the bimodal couplings is degenerate, and, therefore, the process does not have to find a specific unique spin configuration but is entropically attracted to a region with exponentially many ground states and is relaxed once any one out of these many configurations has been reached.

Our results also further reinforce the notion that the relaxation dynamics of SA at these $T=0$ phase transitions is very different from the equilibrium dynamics, where it is known that, with local updates, the exponent governing the ergodic sampling process at fixed finite temperature diverges, $z_{\rm eq}(T)\rightarrow\infty$, when $T\to 0$ liang92 ; wang88 ; katzgraber05 . In contrast, at $T>0$ transitions, in both nonrandom and spin-glass models liu14 ; liu15a ; fernandez16 , the dynamic exponent is finite and takes the same value at equilibrium and in SA analyzed within the KZ hypothesis. Clearly the source of this difference lays in the fact that the equilibrium dynamics is nonergodic in the limit $T\to 0$.

Though the numerical evidence for KZ scaling of the SA dynamics is very strong, we do not have a rigorous theoretical explanation for why it applies, instead of some exponentially slow relaxation dynamics related to naively expected activated scaling. The fact that power-law scaling does hold, in the model studied here as well as in the previously studied case with bimodal couplings rubin17 , must reflect a certain “funnel” structure of the energy landscape where the energy and entropy barriers along the walls down to the global minimum increase sufficiently slowly with the system size. This should be a consequence of the droplet picture in the model with bimodal couplings thomas11 , and also in the case of Gaussian couplings one can construct a similar approximate droplet structure roma10 that may explain the behavior found here.

Given that KZ scaling in the form of data collapse onto a common scaling function is observed, a surprising behavior found here for the Gaussian couplings is that the scaling function for the EA order parameter does not appear to have a power-series expansion for small values of the standard KZ variable $vL^{z+1/\nu}$; instead the data show that the the scaling function has a Taylor expansion in the related variable $Lv^{1/(z+1/\nu)}=L/\xi_{v}$. This indicates a break-down of standard perturbative mechanisms behind KZ scaling in the low-velocity limit, which have been worked out for quantum many-body systems under Hamiltonian dynamics (quantum annealing) degrandi11 and have been shown to be applicable also for stochastic SA dynamics of classical systems liu14 . While the reasons for the non-perturbative behavior found here are not presently clear and deserve further study, one possibility is the proliferation of excited states nearly degenerate with the unique ground state, which may shrink the radius of convergence of the perturbation series to zero in the thermodynamic limit. How these non-perturbative effects lead to analytic behavior in the new scaling argument $L/\xi_{v}$ is not clear and is an important question for further study. Our result for the excess energy are not sufficiently precise to analyze the low-velocity corrections in that case.

Acknowledgements.

We thank David Huse and Anatoli Polkovnikov for helpful discussions. This work was supported by the NSF under Grant No. DMR-1410126 (NX, SJR, and AWS), by MOST in Taiwan through Grants No. 104-2112-M-002-022-356-MY3 and 105-2112-M-002-023-MY3 (KHW and YJK), and by Boston University’s Undergraduate Research Opportunities Program (SJR). Most of the computations were done on the National Center for High-performance Computing’s Formosa 5 Cluster (Taiwan) and some of them were carried out on Boston University’s Shared Computing Cluster.