{'fraction_non_alphanumeric': 0.030834362321397073, 'fraction_numerical': 0.011536426177456342, 'mean_word_length': 4.654498304408538, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Correlations between the dynamics of parallel tempering and the free-energy landscape in spin glasses

We present the results of a large-scale numerical study of the equilibrium three-dimensional Edwards-Anderson Ising spin glass with Gaussian disorder. Using parallel tempering (replica exchange) Monte Carlo we measure various static, as well as dynamical quantities, such as the autocorrelation times and round-trip times for the parallel tempering Monte Carlo method. The correlation between static and dynamic observables for 5000 disorder realizations and up to 1000 spins down to temperatures at 20% of the critical temperature is examined. Our results show that autocorrelation times are directly correlated with the roughness of the free-energy landscape.


I. INTRODUCTION
Systems in statistical physics with rough free-energy landscapes, such as spin glasses, are notoriously difficult to equilibrate using Markov chain Monte Carlo methods. The basic reason for this difficulty is that extremely long times are required to surmount free-energy barriers and, thus, to efficiently and fully explore the minima in the free-energy landscape. The parallel tempering algorithm (also known as "replica exchange Monte Carlo") [1][2][3], as well as other multicanonical techniques [4,5], have partially overcome this problem. However, there is little understanding of the relevant time scales for these methods. The purpose of this work is to shed light on the equilibration time of the parallel tempering algorithm and to understand how it is related to equilibrium properties such as the roughness of the free-energy landscape. To this end, we apply parallel tempering to a large ensemble of disorder realizations of the three-dimensional Edwards-Anderson (EA) Ising spin glass [6,7]. Our main results are that equilibration times vary significantly from one disorder realization to another and that this variation is correlated with the structure of the free-energy landscape. The fact that equilibration times are broadly distributed has been previously noted in a different context (see, for example, Refs. [8] and [9]) but the extent of variation was not fully appreciated because different measures of equilibration time were used in past studies. The correlation of dynamic time scales for parallel tempering and equilibrium properties has not been previously recognized, again because of the use of insensitive measures of dynamic time scales. In this study we compute autocorrelation functions of observables to measure equilibration times whereas previous studies have primarily relied on the round-trip time and related quantities. Although our study is focused on parallel tempering and the Edwards-Anderson Ising spin-glass model, we believe that our qualitative results are likely to be relevant both to other disordered systems with rough free-energy land-scapes, as well as other multicanonical simulation techniques.
Parallel tempering was independently introduced by Geyer [1] and Marinari and Parisi [2], as well as Hukushima and Nemoto [3]. Furthermore, an earlier algorithm with many of the same essential features was introduced by Swendsen and Wang [10]. In parallel tempering many copies of the system at different temperatures are simulated in parallel. Each copy (henceforth referred to as replica) is simulated using a standard, single-temperature Monte Carlo scheme such as the Metropolis or heat bath algorithms [11]. The set of temperatures is chosen so it includes both high temperatures (where rapid equilibration is possible using the single-temperature Monte Carlo method) and low temperatures of interest (where equilibration is not feasible using the standard algorithm). In addition to the singletemperature Monte Carlo moves, there are also replica exchange moves in which replicas at neighboring temperatures swap their temperatures. Replica exchange moves permit replicas to diffuse in temperature space [12]. The diffusion of replicas in temperature space permits an enormous speed-up relative to single-temperature Monte Carlo. Instead of directly surmounting the high free-energy barriers present at low temperatures, barriers are indirectly surmounted via the following process: A replica in one free-energy well at the lowest temperature diffuses to the highest temperature where it easily moves to another well and then diffuses back to the lowest temperature. This process is called a "round trip" and it allows parallel tempering to reduce exponential time scales for surmounting free-energy barriers to power-law time scales [13].
In the context of simple models with rough free-energy landscapes it has recently been shown [13] that the equilibration time of parallel tempering depends strongly on the structure of the free-energy landscape. If there is a single dominant minimum in the landscape and no firstorder transition within the range of temperatures, then the equilibration time is short. On the other hand, the presence of several nearly degenerate minima or the existence of a first-order transition leads to much longer equilibration times. These results in simple model landscapes prompted us to look for similar phenomena in the Edwards-Anderson model where different disorder realizations may have quite different free-energy landscapes. In this study, rather than measuring the free-energy landscape directly, we used the spin overlap (order parameter) distribution as a proxy. The overlap distribution is closely related to the free-energy landscape but is substantially easier to measure in simulations. Based on the results from simple model free-energy landscapes, we hypothesize that when the overlap distribution has a complex structure, implying many free-energy minima the time scale for parallel tempering will tend to be longer than when the overlap distribution is simple. Our results confirm this hypothesis.
This study is part of a larger project whose goal is to understand the nature of the low-temperature phase of finite-dimensional spin glasses. The low-temperature phase of the Edwards-Anderson model is poorly understood and hotly debated. There are several competing theories for the low-temperature behavior of the model [14][15][16][17][18][19][20][21][22][23] and even though many large-scale investigations have been conducted [22][23][24][25][26], a conclusive theory correctly explaining all phenomena has not been agreed on. Our recent simulations [27] point toward a two-pure-state picture such as the droplet picture [14,15,[28][29][30] or the chaotic pairs picture [19,20,31]. One motivation for this work was to verify that the ensemble of disorder realizations used in Ref. [27] was indeed equilibrated. Therefore, we have measured both equilibration times and overlap distributions for all disorder realizations and correlated these quantities. We also demonstrate that the conservative use of the equilibration criterion introduced in Ref. [23] is sufficient to ensure that nearly all disorder samples are equilibrated.
Our results suggest both a method for improving parallel tempering and a warning when using it in spin-glass simulations. On the one hand, we find that many disorder realizations have quite short equilibration times. Thus it might be useful to implement an adaptive scheme where some disorder realization are simulated for shorter times than other realizations. State-of-the-art parallel tempering simulations of equilibrium spin glasses require huge amounts of CPU time because of the difficulty of reaching equilibrium and the need for a large ensemble of disorder realizations. Thus, such an adaptive scheme has the potential for large savings in CPU time. On the other hand, the fact that the overlap distribution (a quantity directly related to the controversy surrounding the low-temperature phase of spin glasses) is significantly correlated with the dynamics of the algorithm serves as a warning that one must be extremely careful to ensure that essentially all samples are well equilibrated in order to avoid systematic errors in measuring disorderedaveraged equilibrium properties.
The paper is organized as follows: In Sec. II we in-troduce the model, the parallel tempering algorithm, as well as the measured observables. In Sec. III we show our results, followed by conclusions.

A. The Edwards-Anderson Ising spin glass
We study the three-dimensional EA Ising spin-glass model, defined by the energy function The sum is over the nearest neighbors on a simple cubic lattice with periodic boundary conditions and side length L. The Ising spins s i take values ±1 and the interactions J ij are quenched random couplings chosen from a Gaussian distribution with zero mean and unit variance.

B. Parallel tempering Monte Carlo
Parallel tempering, also known as replica exchange Monte Carlo, is a powerful tool for simulating systems with rough free-energy landscapes [12,13,32] with applications across many disciplines. To date, it is the most efficient method for simulating spin glasses and other disordered systems in more than two dimensions at low temperatures where dynamics are slow.
The algorithm works as follows: Several copies (hereafter also referred to as replicas) of the system with the same disorder are simulated in parallel. Each replica is simulated at a different temperature using a standard Markov chain Monte Carlo (MCMC) technique such as the Metropolis or heat-bath algorithms [11,33,34]. In this work we use the heat bath algorithm because of its better performance at low temperature. The set of temperatures is chosen to span both the low temperatures of interest, where equilibration is not feasible using single-temperature MCMC methods, and high temperatures where equilibration is fast. In addition to the MCMC sweeps of each temperature, there are also replica exchange moves between replicas. A proposed replica exchange move involves swapping the temperatures of two replicas at neighboring temperatures. A replica exchange move is accepted with probability p, where Here β = 1/T is the inverse temperature of one replica and β ′ the inverse temperature of the neighboring replica, and E and E ′ are the corresponding energies of the two replicas. Equation (2) ensures that the entire Markov chain, including both single-temperature MCMC sweeps and replica exchange moves satisfies detailed balance and converges to equilibrium. Replica exchange allows replicas to diffuse from the lowest temperature to the highest temperature and back again. These round trips enhance mixing and greatly reduce equilibration times for rough energy landscapes. The parallel tempering algorithm has free parameters that include the set of temperatures and the ratio of single-temperature sweeps to replica-exchange sweeps. Optimizing parallel tempering therefore requires an appropriate choice of these parameters. The choice of the number of replicas N T and their temperatures involves the following trade-off: If there are many closely-spaced temperatures, the energy distribution between adjacent temperatures overlaps strongly and the probability that a proposed replica exchange is accepted, see Eq. (2), is large. This suggests that the more replicas one uses, the better, apart from the extra computational work involved in carrying out single-temperature sweeps of each replica. However, the motion of replicas in temperature space is diffusive such that the time scale for a round-trip scales approximately as the square of the number of temperatures. At least for simple model systems, parallel tempering is optimized if the number of temperatures scales as the square root of the number of spins [13]. The average time to complete one round trip is often used to characterize the performance of parallel tempering. Choosing a set of temperatures that minimizes the round trip time is one of the ways proposed to optimize parallel tempering [12,[35][36][37].

C. Simulation parameters
For the parallel tempering algorithm we use N T = 16 temperatures between T = 0.2 to 2.0 [38]. This set of temperatures was chosen heuristically in Refs. [23] and [24] to perform well for L = 6 and 8. The corresponding average acceptance fractions for replica exchange moves are plotted in Fig. 1. It is likely that the simulations would be more efficient if more replicas were used near and above the critical region (T c ≈ 0.96 [25]) where the acceptance fractions are small. However, we believe that as long as the swapping probability is nonzero, our results will not change qualitatively. A bottleneck merely slows down the diffusion process but does not prevent it. For L = 10, the small acceptance fractions suggest that a larger number of temperatures might be more efficient. However, our goal is to understand correlations between static and dynamic quantities in parallel tempering for a given system size and so we have not sought to optimize the algorithm for each system size. We believe that the qualitative results obtained using our parallel tempering parameters would also hold for other reasonable values of the parameters.
By placing the hottest system temperature at T ≈ 2T c , we ensure that equilibration happens quickly, and we use one exchange sweep for each heat-bath sweep (N = L 3 attempted spin updates). To calculate the spin overlap [Eq. (3)] we use two copies of the system at each temperature; therefore, for each sample, we simulate two sets of replicas independently from each other.
As stated before, we employ one replica exchange sweep after one heat-bath sweep. A heat-bath sweep corresponds to sequentially attempting to update each spin at each replica once using the heat-bath algorithm. An exchange sweep corresponds to randomly choosing N T −1 pairs of adjacent replicas and proposing exchanges for each. Both types of sweeps together make up one parallel tempering sweep. When no confusion arises, we will call this unit of time simply a sweep.
For each system size we simulate approximately 5000 disorder realizations or samples. The simulations are equilibrated for at least 2 24 sweeps for L = 6 and 2 27 sweeps for L = 8 and 10, and then measurements are done for the same number of sweeps; see Table I for details. A fraction of samples required longer runs to meet the equilibration criteria discussed in Sec. II D 2. For samples with very short equilibration times, we performed runs to obtain fine-grained autocorrelation information; see Sec. II D 2. In total, the data collection took approximately 140 CPU years on 12-core AMD Opteron 6174 CPUs.

Overlap distributions
The focus of the paper is to relate static equilibrium properties of individual spin-glass samples to the dynamics of the parallel tempering algorithm acting on that sample. The primary quantity that we measure to study both equilibrium properties and the dynamics of parallel tempering is the spin overlap q, where N = L 3 is the number of spins and the superscripts (1) and (2) indicate two independent copies of the system with the same disorder. The thermal and disorder average of the overlap is an order parameter for the system and the probability distribution of the thermally averaged overlap P J (q) for a given sample J reveals aspects of the free-energy landscape for that particular sample. The behavior of P J (q) for large systems and at low temperatures is one of the major open questions in the theory of spin glasses. Overlap distributions vary widely from sample to sample, as can be seen in the three examples shown in Fig. 2. Although there is no direct mapping between the free-energy landscape and P J (q), it is clear that numerous peaks in P J (q) imply numerous minima in the free-energy landscape. For example, samples corresponding to Figs. 2(b) and 2(c) have a more complicated free-energy landscapes than Fig. 2(a). Our central hypothesis is that samples with a more complicated free-energy landscape tend to have longer dynamic time scales. From P J (q) we define I J (q 0 ) for different q 0 , where I J (q 0 ) is the weight of P J (q) in an interval near the origin and it serves as an approximate measure of the complexity of P J (q). For example, I J (0.2) ≈ 0 for the samples shown in Figs. 2(a) and 2(b), while I J (0.2) is large for the sample shown in Fig. 2(c). We compute I J (q 0 ) for eight values of q 0 from q 0 = 0.2 to 0.9. It is worth noting that I J (q 0 )-especially with q 0 = 0.2-has been extensively used in studies of the low-temperature equilibrium properties of the EA model [23,26].

Characteristic time scales
We measure two time scales for parallel tempering from the autocorrelation function of the overlap at the lowest temperatures. For an arbitrary observable A, the autocorrelation function Γ A (t) is defined via Here, A(t) is the observable measured at time t and · · · represents a thermal average [39]. Integrated and exponential autocorrelation times are computed from the autocorrelation function. The integrated autocorrelation time τ A int is the integral of the autocorrelation function or, for discrete time, the sum of the autocorrelation function, The integrated autocorrelation time is the time needed for two subsequent measurements of A to decorrelate. It is a lower bound on the time needed to reach equilibrium from an arbitrary initial condition. To calculate τ A int care must be taken when truncating the sum in Eq. (6) to avoid large statistical errors (explained below).  We obtain the noise floor for each run length by visual inspection of Γ |q| (t) for many samples. We note that this approach introduces some slight systematic bias in our estimates of the integrated autocorrelation times. However, the main emphasis of this paper is to understand the very large sample-to-sample variations in time scales and, compared to these variations, the errors introduced by our truncation protocol are small. We have computed autocorrelation functions for both |q| and q and from these we have obtained τ |q| int and τ q int , respectively. To do this, we first compute a Fourier transform of the data and then invert the function to obtain the autocorrelation function. This method is considerably faster than a direct calculation. For L = 6, the simulations are run for 2 24 sweeps, recording data for the computation of the autocorrelation functions every 10 sweeps. The noise floor for these runs is set to Γ = 0.005. Some samples have τ |q| int values of order 10 or less. For these, shorter data collection runs of 2 18 sweeps have been done starting from an equilibrated spin con-figuration stored at the end of the longer run. For these short runs data are collected every sweep to accurately measure the autocorrelation function up to time 10 and the full autocorrelation function is patched together from the short and long runs. The noise floor for these shorter runs is set to Γ = 0.01. If the noise floor is reached before t = 10, only the autocorrelation function generated by the short run is used. For both L = 8 and L = 10 the simulations are run for at least 2 27 sweeps, recording data every 100 sweeps with shorter runs also needed for samples with small τ |q| int values. The noise floor used is the same as for L = 6.
The autocorrelation function is always calculated up to a maximum time lag of 1% of the total run length. Therefore, for some samples and, in particular for L = 10, the noise floor was not reached by this time and longer runs were necessary to obtain good statistics. For these L = 10 samples we simulate up to 2 30 sweeps. 32 samples (approximately 0.6%) stayed above the noise floor even for the longest runs. To prevent a biasing of the results, these were not included in the analysis. However, it would be interesting to study these samples in detail in the future.
There is another method for truncating the sum in Eq. (6) introduced by Madras and Sokal [40] where the upper limit of the sum is determined recursively. An initial upper limit is chosen and τ A int is computed for that value, and the new upper limit is determined as some factor (6 is suggested in Ref. [40]) times the current τ A int . The estimated value of τ A int is obtained when this process has converged. We experimented with this approach but did not find that it converged in a consistent way across all samples with very different autocorrelation functions. However, we note that our approach should have systematic errors no larger than those resulting from the method of Madras and Sokal with the upper limit set to 6τ The second time scale that we consider is the exponential autocorrelation time τ A exp , defined by the asymptotic exponential decay of the autocorrelation function, Except for observables orthogonal to the slowest mode of the Markov chain, τ A exp is the characteristic time of the slowest mode and is thus the exponential time scale to reach equilibrium from an arbitrary initial state [41].
In principle, τ A exp is the most important time scale for studying disordered systems because the main difficulty is reaching thermal equilibrium (controlled by τ A exp ) rather than collecting enough uncorrelated data at equilibrium (controlled by τ A int ). Unfortunately, calculating τ A exp for a large number of samples is difficult because it requires an automatized fitting process. For some samples, such as the one shown in Fig. 3(a), the autocorrelation function decays nearly exponentially over most of the observable time range, allowing for a precise fitting and extraction of τ |q| int . However, for many samples the autocorrelation function is not exponential for the measured times, as can be seen in Fig. 3(b) and, for these samples, an automatic fitting procedure would not be reliable. For that reason, the time scale that we correlate with the static properties of the EA spin glass is the integrated autocorrelation time. We have studied τ |q| exp by hand for a small subset of the data set and find that for those samples where a reasonable fit is possible, τ |q| exp is in the range of 1 to 15 times τ |q| int . Most, but not all, samples share the behavior seen in the examples in Figs. 3 that there is an initial sharp decline before the (approximately) exponential decay sets in. This sharp decline explains why τ |q| exp is typically larger than τ |q| int . Finally, we also measured the round-trip time τ RT for each sample, defined as the average time after equilibration is achieved, for a replica to diffuse from the lowest to the highest temperature and then back to the lowest temperature again [12]. Figure 4 shows the probability distribution for τ |q| int for L = 8. In agreement with previous studies [8,9,25], the figure reveals that the equilibration times for spin glasses are very broadly distributed. The distributions for L = 6 and 10 are similar although with different logarithmic means and standard deviations, as shown in Tables II  and III, respectively. Tables II and III also show the logarithmic mean and standard deviation of τ q int . The values for τ q int are similar to those for τ |q| int , however, the standard deviations are somewhat smaller. Figure 5     tained for the other system sizes studied). This suggest that when the free-energy landscape is complex, timescales are large, as one would naively expect.

III. RESULTS
To quantify the correlations between τ |q| int and I J (q 0 ), we use the Pearson correlation coefficient r [42], whose estimator for a finite data set {(x i , y i )} is where M is the size of the data set and the overbar indicates the average over the data set. If r = 1 there is an exact linear relationship between x and y while a small value of r indicates the absence of any linear dependence.
Pearson r values for log 10 (τ |q| int ) vs log 10 [I J (q 0 )] for our data are shown in Table IV and Fig. 7 for the system sizes studied and different values of q 0 . These results clearly demonstrate that longer parallel tempering time scales are correlated with more complicated overlap distributions. The strength of the correlation is slightly weaker for larger systems. The fact that r increases with q 0 suggests that the presence of additional peaks in P J (q) increases τ |q| int independent of whether those peaks are near the origin.
On the other hand, we find that the round trip time is     (17) 0.126 (20) not correlated with τ |q| int or, by extension, I J (q 0 ). Figure  8 is a scatter plot of τ RT vs τ |q| int for L = 8 (compare with Fig. 5). Each point represents one sample and the diagonal line is τ RT = τ |q| int . It is clear from this figure that τ RT has a much smaller variance than τ |q| int and that the correlation between the two is minimal. The last row of Table IV shows the Pearson r value between log 10 (τ RT ) and log 10 (τ |q| int ) and quantifies the lack of correlation between τ RT and τ |q| int . Note that the authors of Ref. [26] consider a time scale related to the motion of replicas in parallel tempering and also find no correlation with I J (q 0 ).
If the highest temperature is chosen such that its mixing time is one heat bath sweep, then τ RT ≥ τ |q| int because there will be no memory of the spin state of the replica between successive visits to the lowest temperature, as can be seen in Tables II and III, as well as Fig. 8. Figure  10  8 shows that τ RT is much larger than τ |q| int for most samples. One explanation for this is that (for most samples) the spin state of the coldest replica can be decorrelated by diffusion from the lowest temperature to an intermediate temperature and back; full round trips are not needed for decorrelation.
Note that if we repeat the study for τ q int instead of τ |q| int , the correlation values would show the same qualitative trend as for τ |q| int , but they would, on average, be slightly smaller.
Finally, we also test equilibration with the method introduced by Katzgraber et al. in Ref. [23] for short-range spin glasses with Gaussian disorder between the spins. In this method, the average energy per spin where [· · · ] av represents a disorder average and, as before, · · · denotes the Monte Carlo average for a given set of bonds, can be related to the link overlap via an integration by parts over the interactions J ij between the spins. In Eq. (10) N b = dN with d = 3 the space dimension represents the number of bonds in the system. One obtains: where T is the temperature.  shows representative results for L = 8 and T = 0.2, the lowest temperature simulated. At approximately 10 5 Monte Carlo sweeps both data computed directly from the energy per spin (blue squares) and data computed from the link overlap (red circles) agree within error bars. Note that 99% of the τ |q| int values for this system size are smaller than the point marked with the solid vertical (green) line in Fig. 9. Thus, for the simulated time of over 10 8 Monte Carlos sweeps, the data are well equilibrated. This suggests that a conservative use of the equilibration criterion developed in Ref. [23] will guarantee that almost all of the samples are equilibrated.

IV. CONCLUSIONS
We have found one factor that explains the broad distribution of time scales in parallel tempering for spin glasses: The roughness of the free-energy landscape of each individual sample directly affects its equilibration time scales. In accordance with previous observations on the parallel tempering method [13], if there are many minima in the free-energy landscape, a sample requires more time to equilibrate. The equilibrium distribution of the spin overlap serves as proxy for the complexity of the free-energy landscape. Our results show that individual samples need to be tested individually to ensure proper equilibration, as previously suggested in Ref. [9]. However, we show that the equilibration test developed in Ref. [23] is a viable alternative for system sizes currently accessible in simulations if applied conservatively. Furthermore, the parallel tempering round-trip times seem to be unaffected by the complexity of the free-energy landscape and should therefore not be used as a mea-sure of how well a system is equilibrated.
The clear correlations between algorithmic (nonequilibrium) time scales and the equilibrium complexity of the free-energy landscape opens the door for alternate studies of the nature of the spin-glass state, a problem that remains controversial.