## Evidence against a mean-field description of short-range spin glasses revealed through thermal boundary conditions

Wenlong Wang ∗

Department of Physics, University of Massachusetts, Amherst, Massachusetts 01003, USA

## Jonathan Machta †

Department of Physics, University of Massachusetts, Amherst, Massachusetts 01003, USA and Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, New Mexico 87501, USA

## Helmut G. Katzgraber

Department of Physics and Astronomy, Texas A&amp;M University, College Station, Texas 77843-4242, USA Materials Science and Engineering Program, Texas A&amp;M University, College Station, Texas 77843, USA and Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, New Mexico 87501, USA

A theoretical description of the low-temperature phase of short-range spin glasses has remained elusive for decades. In particular, it is unclear if theories that assert a single pair of pure states, or theories that are based on infinitely many pure states-such as replica symmetry breaking-best describe realistic short-range systems. To resolve this controversy, the three-dimensional EdwardsAnderson Ising spin glass in thermal boundary conditions is studied numerically using population annealing Monte Carlo. In thermal boundary conditions all eight combinations of periodic vs antiperiodic boundary conditions in the three spatial directions appear in the ensemble with their respective Boltzmann weights, thus minimizing finite-size corrections due to domain walls. From the relative weighting of the eight boundary conditions for each disorder instance a sample stiffness is defined, and its typical value is shown to grow with system size according to a stiffness exponent. An extrapolation to the large-system-size limit is in agreement with a description that supports the droplet picture and other theories that assert a single pair of pure states. The results are, however, incompatible with the mean-field replica symmetry breaking picture, thus highlighting the need to go beyond mean-field descriptions to accurately describe short-range spin-glass systems.

## I. INTRODUCTION

A plethora of problems across disciplines and, in particular, a wide variety of optimization problems map onto spin-glass-like Hamiltonians [1-6]. As such, despite the fact that only a selected class of disordered magnets, such as LiHoF 4 or Au Fe x 1 -x show this intriguing state of matter, spin glasses have been of great importance across multiple fields including condensed matter physics, evolutionary biology, neuroscience and computer science. Most recently spin glasses have played a pivotal role in the development of new computing prototypes based on quantum bits in both a theoretical, as well as device-centered role. For example, the stability of topologically protected quantum computing proposals [7-9] against different error sources-recently implemented experimentally [10]-heavily relies on spin-glass physics [11-13]. Similarly, the native benchmark problem currently used to gain a deeper understanding of state-of-the-art quantum annealing machines is based on a spin-glass Hamiltonian [14-17]. Given this recent renaissance of spin glasses to benchmark novel algorithms, as well as to develop cutting-edge computing paradigms, it is unsettling that no consensus exists as to whether

∗ Electronic address: wenlong@physics.umass.edu

† Electronic address: machta@physics.umass.edu

mean-field theory, also known as replica symmetry breaking theory (RSB) [2, 3, 18-21], accurately describes the low-temperature phase of these systems.

Spin-glass models are well understood in the meanfield regime where infinite-range interactions dominate. However, in finite space dimensions spin glasses are still poorly understood and have been the subject of a long-standing controversy. In this paper we seek to help resolve this controversy. Using numerical methods, we study the low-temperature phase of the threedimensional Edwards-Anderson (EA) spin glass [22]. In so doing, we introduce two methods that promise to be useful in the study of other disordered systems-thermal boundary conditions and sample stiffness extrapolation.

The controversy concerning the EA model is between two competing classes of theories as to the nature of the low-temperature phase. One proposal, championed by Parisi and collaborators [2, 3, 18-21, 23], is that finitedimensional EA spin glasses behave like the mean-field Ising spin glass, known as the Sherrington-Kirkpatrick (SK) model [24]. Parisi analytically studied the SK model [18-20] and found that the low-temperature phase is characterized by an unusual form of symmetry breaking called replica symmetry breaking (RSB). This RSB solution of the SK model predicts that there is an infinity of pure thermodynamic states and that the overlap distribution of these pure states is not self-averaging. Many features of Parisi's RSB solution of the SK model have now been verified by rigorous mathematical methods [25].

The mean-field or RSB picture for finite-dimensional EA spin glasses asserts that the qualitative features of the SK model also hold for finite-dimensional models so that, in particular, there are infinitely many pure states in the thermodynamic limit.

In contrast to the RSB picture, the main competing class of theories for the three-dimensional (3D) EA model assume that the low-temperature phase consists simply of a single pair of pure states related by the spin-reversal symmetry of the Hamiltonian. The earliest and most widely accepted of these theories is the 'droplet picture' developed by McMillan [26], Bray and Moore [27], and Fisher and Huse [28-30]. The droplet picture asserts that the low lying excitations of the pure states are compact droplets with energies that scale as a power of the size of the droplet. By contrast, the low lying excitations in the RSB picture are space filling objects.

Several features of the original RSB picture for finitedimensional EA models have been mathematically ruled out in a series of papers by Newman and Stein [31-35]. These authors provide two alternative theories for finitedimensional EA models, both of which have infinitely many pure states. The first is a nonstandard RSB picture, similar to the original RSB picture but with a selfaveraging thermodynamic limit. Newman and Stein give heuristic arguments against the nonstandard RSB picture but do not rule it out. On the other hand, the nonstandard RSB picture is promoted as a viable theory for finite-dimensional EA models in Ref. [36]. The second is the 'chaotic pairs' picture. Here there are infinitely many pure states but they are organized in such a way that in each finite volume only a single pair of states related by a global spin flip is seen.

In the following we refer to all pictures that display a single pair of pure states in each large finite volume as two-state pictures. Therefore, the droplet and chaotic pairs pictures are both two-state pictures within this definition. Note that for the droplet model it is the same pair of states in every volume while for chaotic pairs a different pair of states is manifest in each volume.

Parallel to these analytical efforts, many computational studies have been aimed at distinguishing between the two classes of theories (see, for example, Refs. [3745]). Unfortunately, computational methods have been difficult to apply to spin glasses. The fundamental questions concern the limit of large system sizes, however, attempts to extrapolate to large sizes have not been conclusive because the range of sizes accessible to simulations at low temperatures is quite small and, for fixed size, the variance between samples for many observables is quite large. Thus, a straightforward extrapolation to large sizes based on mean values of observables can be misleading. Computational studies have yielded a confusing mixture of results: Some point to the RSB picture, some to a two-state picture, and some to a mixed scenario, known as the the 'trivial nontrivial' (TNT) picture described in Refs. [38, 42, 46]. Recently, there have been efforts to analyze statistics other than simple disorder averages [44, 47-51] but these methods have not been definitive either and so the controversy continues.

Here we introduce two related innovations to more effectively extrapolate from small system sizes to the large system-size limit. First, we employ thermal boundary conditions instead of the usual periodic boundary conditions. In d space dimensions, thermal boundary conditions allow the 2 d combinations of periodic or antiperiodic boundary conditions each to appear with the correct Boltzmann weights. The idea is to let the system choose boundary conditions that minimize the presence of domain walls and thus finite-size effects. Second, based on thermal boundary conditions, we define a spin stiffness measure for each sample. We show, as expected for a low-temperature phase, that the sample stiffness becomes large as the system size becomes large. We then study the behavior of the system in the limit of large sample stiffness and relate the system's behavior for large sample stiffness to the large system-size limit. During this process we also obtain new measurements of the spin stiffness exponent at nonzero temperatures. These generic techniques promise to be of broad utility in understanding disordered systems in statistical mechanics, not just the EA spin-glass model. Furthermore, using this approach we conclude that a two-state picture best describes the low-temperature phase of the threedimensional Edwards-Anderson Ising spin glass.

The paper is structured as follows. In Sec. II we introduce the studied model, as well as thermal boundary conditions. Section III describes the implementation of population annealing Monte Carlo used in this study, followed by the measured quantities in Sec. IV. Results are presented in Sec. V, followed by a discussion in Sec. VI and concluding remarks.

## II. THE EDWARDS-ANDERSON MODEL IN THERMAL BOUNDARY CONDITIONS

## A. Edwards-Anderson model

We study the three-dimensional Edwards-Anderson Ising spin-glass model [22]. The model is defined by the Hamiltonian

$$\mathcal { H } = - \sum _ { \langle i, j \rangle } J _ { i j } S _ { i } S _ { j },$$

where S i ∈ {± } 1 are Ising spins and the sum is over nearest neighbors on a cubic lattice of linear size L . The random couplings J ij are chosen from a Gaussian distribution with zero mean and unit variance. A set of couplings J = { J ij } defines a disorder realization or 'sample.'

## B. Thermal boundary conditions

Most Monte Carlo simulations of spin systems are performed with periodic boundary conditions (PBC) be-

cause it is often assumed that periodic boundary conditions yield the mildest finite-size correction. Free boundary conditions are sometimes also employed [39] as are antiperiodic boundary conditions in one direction for the purpose of measuring spin stiffness. In this work we consider thermal boundary conditions (TBC). Thermal boundary conditions include the set of all 2 d choices of periodic or antiperiodic boundary conditions in d spatial dimensions. This means that for three space dimensions ( d = 3) we have 8 possible choices. For example, one of the 8 elements in the set of boundary conditions is 'periodic in the x -direction and antiperiodic in the y and z directions.' For each boundary condition ζ , there is a free energy F ζ and the probability distribution for spin states in TBC is the weighted mixture of the eight boundary conditions with weights e -βF ζ . An equivalent way to describe thermal boundary conditions is to say that the eight boundary conditions are annealed so that each spin configuration together with a boundary condition appears with its proper Boltzmann weight.

The motivation for using thermal boundary conditions for spin glasses can be explained by considering two simpler examples-the ferromagnetic and antiferromagnetic Ising models on a square lattice with lattice size L an odd number. For the ferromagnetic Ising model periodic boundary conditions in all directions are natural and appropriate because they do not induce domain walls in the ordered phase. However, for the antiferromagnetic Ising model, periodic boundary conditions will induce domain walls and the observables for finite systems will have strong finite-size corrections. The natural boundary conditions for the antiferromagnet with L odd are antiperiodic in all directions. Now, suppose we are asked to simulate an Ising model but we are not told whether it is a ferromagnet or antiferromagnet. If we use thermal boundary conditions then we will automatically choose the natural boundary conditions independent of which model we have been given, namely periodic in all directions if the system is a ferromagnet and antiperiodic in all directions if the system is an antiferromagnet. The other boundary conditions will induce domain walls and therefore have higher free energies. The difference in free energy between thermal boundary conditions and any of the domain-wall-inducing boundary conditions scales as L θ where θ is the spin stiffness exponent. For the Ising model (either ferromagnetic or antiferromagnetic) in the low-temperature phase, θ = d -1 ≥ 0, and, even for modest system sizes, thermal boundary conditions are essentially the same as the single natural boundary condition because all unfavorable choices are suppressed.

While one can a priori determine the optimal boundary conditions for simple systems such as ferromagnets and antiferromagnets, the same is not true for spin glasses. For a given sample, a single boundary condition such as PBC may induce domain walls and induce large finite-size effects. The motivation for using thermal boundary conditions is thus the same as for the simple (anti)ferromagnetic example discussed above. Because we do not know which of the eight periodic/antiperiodic boundary conditions fits the sample best, we simply let the system choose by minimizing the free energy.

At zero temperature, thermal boundary conditions correspond to selecting from among the 2 d boundary conditions those with the lowest energy ground states. These boundary conditions have been employed with exact algorithms for finding ground states of two-dimensional spin glasses [52, 53]. Thomas and Middleton [53] call thermal boundary conditions 'extended' boundary conditions and argue, as we do, that these boundary conditions minimize finite-size effects. Similar ideas but using periodic and antiperiodic boundary conditions in a single direction are discussed in [54-57].

For the mathematical statistical physicist the difficulties produced by spurious domain walls are avoided by using 'windows:' Consider a very large system of linear size L and a large window inside this system of linear size glyph[lscript] such that 1 glyph[lessmuch] glyph[lscript] glyph[lessmuch] L . Then any domain walls induced by the 'bad' boundary conditions will almost surely lie outside the window and observables measured within the window will not be influenced by the domain wall. By collecting data only inside the window the results are then independent of the boundary conditions. Unfortunately, the computational statistical physicist does not have the luxury of collecting equilibrated data in this way for spin glasses where attainable system sizes deep within the low-temperature phase do not exceed, for example, L ≈ 10 in three space dimensions.

Thermal boundary conditions may also be used to measure the spin stiffness exponent θ by comparing the free energy of TBC with other boundary conditions. For example, for spin glasses it is sufficient to compare the free energy for TBC with that for PBC. This approach is expected to yield the same exponent but a different prefactor for the spin stiffness as compared to the standard method of taking the absolute value of the free energy difference between periodic and antiperiodic boundary conditions.

## III. METHODS

We use population annealing Monte Carlo [58-60] to simulate the EA model. The most common method for large-scale spin-glass simulations is parallel tempering Monte Carlo [61], however, parallel tempering Monte Carlo is not well suited to thermal boundary conditions and is not easily used to measure free energies. In contrast, population annealing is able to simulate TBC and accurately measure free energies in a straightforward way.

Population annealing is related to simulated annealing where the temperature of the system is lowered in a step-wise fashion following a predefined annealing schedule. In each annealing step a Markov chain Monte Carlo algorithm is applied to the system at the current temperature in the schedule. In population annealing a large population of replicas of the system are simultaneously

annealed from high to low temperature. However, simulated annealing is designed for finding only ground states and it does not correctly sample equilibrium states at the temperatures traversed during the annealing schedule. Population annealing corrects this deficiency by adding a resampling step that ensures that the population of replicas at every temperature is a Gibbs ensemble at that temperature. Population annealing is an example of a sequential Monte Carlo algorithm [62] and it converges to the correct equilibrium distribution as the population size increases. It is thus well suited to parallel computation and our implementation uses OpenMP.

The algorithm works as follows. Let R 0 be the initial size of the population of replicas of the system. In our implementation of population annealing, each replica is initialized independently at infinite temperature T ( β = 1 /T = 0). Each replica has the same set of couplings. For thermal boundary conditions, 1 / 8 of the replicas are assigned to each of the 8 boundary conditions. The temperature of the population is now cooled from β = 0 in a sequence of steps to a target temperature β 0 . The annealing step from β to β ′ ( β ′ &gt; β ) consists of two stages: The first stage is resampling and the second stage is the application of the Metropolis algorithm at inverse temperature β ′ . In the resampling step, some replicas are eliminated and others are duplicated. For TBC, when a replica is copied, its boundary condition is copied with it.

The resampling step works as follows. Suppose we have ˜ R β replicas that represent an equilibrium ensemble at inverse temperature β and we want to lower the temperature to β ′ &gt; β . The ratio of the statistical weight at β to β ′ for replica j , with energy E j is exp[ -( β ′ -β E ) j ]. In principle this factor should represent how many copies to make of the system. However, this ratio is typically larger than unity. In order to keep the population size roughly fixed we need to normalize the ratio. First compute normalized weights τ j ( β, β ′ ) whose sum over the ensemble is R 0 ,

$$\tau _ { j } ( \beta, \beta ^ { \prime } ) = \left ( \frac { R _ { 0 } } { \tilde { R } _ { \beta } } \right ) \frac { \exp { [ - ( \beta ^ { \prime } - \beta ) E _ { j } ] } } { Q ( \beta, \beta ^ { \prime } ) }, \quad ( 2 ) \quad \text{when} \quad \text{that} \quad \text{$\quad$}$$

where Q is the normalization given by

$$Q ( \beta, \beta ^ { \prime } ) = \frac { \sum _ { j = 1 } ^ { \tilde { R } _ { \beta } } \exp \left [ - ( \beta ^ { \prime } - \beta ) E _ { j } \right ] } { \tilde { R } _ { \beta } }. \quad \ \ ( 3 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

The new population at temperature β ′ is obtained by differential reproduction. The number of copies n j of replica j is either the floor (greatest integer less than) glyph[floorleft] τ j glyph[floorright] or ceiling glyph[ceilingleft] τ j glyph[ceilingright] with probabilities glyph[ceilingleft] τ j glyph[ceilingright] -τ j and τ j -glyph[floorleft] τ j glyph[floorright] , respectively. Note that this choice ensures that the expectation of n j is τ j . It also minimizes the variance of n j among all integer probability distributions with this expectation. (Note that it is possible to have n j = 0 for replica j .) In our implementation, the size of the population at each temperature ˜ R β is variable but stays close to the target

TABLE I: Parameters of the numerical simulations for different system sizes L and periodic (PBC), as well as thermal (TBC) boundary conditions. R 0 represents the number of replicas, 1 /β 0 is the lowest temperature simulated, N T the number of temperatures used in the annealing schedule, N S the number of sweeps per temperature, and M the number of samples. M PBC is the number of hard samples for periodic boundary conditions and M TBC the number of hard samples for thermal boundary conditions.

|   L | R 0 PBC   | R 0 TBC   | 1 /β 0   |   N T |   N S |    M |   M PBC |   M TBC |
|-----|-----------|-----------|----------|-------|-------|------|---------|---------|
|   4 | 5 10 4    | 5 10 4    | 0 . 200  |   101 |    10 | 4941 |       0 |       0 |
|   6 | 2 10 5    | 2 10 5    | 0 . 200  |   101 |    10 | 4959 |       0 |       0 |
|   8 | 5 10 5    | 5 10 5    | 0 . 200  |   201 |    10 | 5099 |       5 |      33 |
|  10 | 10 6      | 2 10 6    | 0 . 200  |   301 |    10 | 4945 |     286 |     291 |
|  12 | 10 6      | 3 10 6    | 0 . 333  |   301 |    10 | 5000 |     533 |     386 |

value R 0 . In the next stage of the annealing step, every member of the new population is subject to N S sweeps of the Metropolis algorithm at inverse temperature β ′ . Our implementation of population annealing follows a schedule of N T inverse temperatures β that are evenly spaced between β = 0 and β 0 . The system sizes and temperatures in the simulations, together with the parameters of the population annealing simulations are shown in Table I. The definition of hard samples is given below. Note that although the number of sweeps per temperature N S is small, the total number of sweeps, R N N 0 S T , is large and comparable to the number of sweeps performed in parallel tempering simulations. In population annealing, equilibration results from large values of R 0 and is guaranteed in the limit R 0 →∞ for fixed N S and N T .

The free energy can be estimated [59] from the normalization factors Q β,β ( ′ ) defined in Eq. (2) according to

$$- \, \beta _ { k } \tilde { F } ( \beta _ { k } ) = \sum _ { \ell = N _ { T } - 1 } ^ { k + 1 } \ln Q ( \beta _ { \ell }, \beta _ { \ell - 1 } ) + \ln \Omega, \quad ( 4 )$$

where Ω is the number of configurations of the system so that for N Ising spins, Ω = 2 N for PBC and Ω = 2 ( N + ) d for TBC in d dimensions.

Because population annealing has not been used before for large-scale simulations in statistical physics, we did a careful comparison to data previously obtained using parallel tempering [44, 63]. We measured observables for the same set of samples studied in Refs. [44, 63], comprising approximately 5000 samples for each system size L with periodic boundary conditions. We found no statistical difference between population annealing and parallel tempering for any disorder averaged observable (see Sec. V B for a comparison of one observable). A detailed comparison between population annealing and parallel tempering will be presented in a subsequent publication [64].

The convergence to equilibrium of population annealing for each sample can be quantified using the family

entropy . Define family i as the set of replicas at some low temperature that are descended from replica i at the highest temperature. In practice most families are empty sets. Let η i be the fraction of the population in family i , i.e., the fraction of the population at the low temperature that is descended from replica i in the initial, high-temperature population. Then the family entropy S f is given by

$$S _ { f } = - \sum _ { i } \eta _ { i } \log \eta _ { i }. \quad \ \ ( 5 ) \stackrel { \dashv \dashv } { \text{ref} } \\ \stackrel { \text{sta} } { \text{sct} }$$

The exponential of S f is an effective number of families. For example, if there are k surviving families all of the same size then η i = 1 /k for each surviving family and e S f = k . In practice, the family sizes are exponentially distributed.

Since each family has an independent history during the simulation, e -S / f 2 is a conservative measure of statistical errors. As discussed in Ref. [64], e -S f is a reasonable measure of systematic errors. In our simulations, we require that e -S f &lt; 0 01 for every disorder sample. . For hard samples that do not meet this requirement with the standard population size given in Table I we increased the population size until this equilibration criterion was met. The numbers of hard samples M PBC and M TBC for periodic and thermal boundary conditions, respectively, are given in Table I. Most hard samples were equilibrated using 5 runs with R 0 = 3 × 10 , which were then com6 bined using weighted averaging [59] for a total population of 1 5 . × 10 . 7 For the hardest samples of size L = 10 and 12, population sizes up to 10 8 were required to meet the equilibration criterion.

## IV. MEASURED QUANTITIES

## A. Free energy, ground-state energy and spin stiffness

Using Eq. (4) we measure the free energies, F TBC J and F PBC J , for each sample J in thermal (TBC) and periodic (PBC) boundary conditions, respectively. We also measure the ground-state energies E TBC J and E PBC J for each sample in both TBC and PBC, respectively. We compute the ground-state energy by taking the minimum energy in the population at the lowest temperature ( T = 0 2 . glyph[lessmuch] T c ). We report on a careful study of the ground-state calculation in a separate paper [65]. There we show that the average ground-state energy agrees with other methods, that multiple runs always yield the same ground state, and that for a small number of the hardest samples we find agreement with exact branch and bound methods.

The traditional measure of spin stiffness is the difference between the free energy, or at zero temperature, the ground-state energy of two different boundary conditions-usually periodic and antiperiodic in a single direction with periodic boundary conditions in all other directions. For spin glasses, this quantity may be of either sign and the absolute value must be taken before performing the disorder average. Here we consider the free energy (ground-state energy) difference between thermal boundary conditions and periodic boundary conditions. This quantity is nonnegative because periodic boundary conditions are contained in the TBC ensemble of boundary conditions so no absolute value needs to be taken. We refer to ∆ F as the disorder average free energy (groundstate energy) difference between TBC and PBC. The scaling of ∆ F with system size L defines the stiffness exponent θ ,

$$\Delta F \sim L ^ { \theta }.$$

We measure θ at T = 0, 0 2, and 0 42 by fitting to this . . equation.

The free energy for each boundary condition in the TBC ensemble can be measured using an analog of Eq. (4) by partitioning Q into its 8 boundary condition components but we did not collect data to do this measurement. Instead, we estimate the ratio of the free energy of the dominant boundary condition to the free energy of all the other boundary conditions combined. Let f J be the fraction of the population in the boundary condition with the largest population in sample J . The quantity λ J ,

$$\lambda _ { \mathcal { J } } = \log \frac { f _ { \mathcal { J } } } { ( 1 - f _ { \mathcal { J } } ) },$$

is an estimator of the free-energy difference (times -β ) between the dominant boundary condition and all other boundary conditions in sample J . Note that λ J is a measure of the stiffness of sample J . If only one boundary condition dominates the ensemble of boundary conditions it means that inserting a domain wall is very costly and the sample is stiff while if λ J is small, the domain walls induced by changing boundary conditions have little cost and the sample is not stiff. Note that, in principle, all boundary conditions could have equal weight so λ ≥ -log 7.

## B. Order parameter distribution

The order parameter for spin glasses is obtained from the spin overlap, q defined by

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

where the superscripts '(1)' and '(2)' indicate two statistically independent spin configurations chosen from the Gibbs distribution. Let P J ( q ) be the overlap distribution for sample J and let P q ( ) be the disorder average of the overlap distribution. In population annealing, the pairs of independent spin configurations used in Eq. (8)

are chosen randomly from the population of replicas with the restriction that the two replicas are from different families. This ensures that the spin configurations are independent.

Two-state pictures make very different predictions from the RSB picture for P J ( q ) for the low-temperature phase of the EA model in the infinite-volume limit. If there is a single pair of pure states then P J ( q ) consists of two δ functions at ± q EA and, of course, P q ( ) after disorder averaging is the same. Here q EA is the EdwardsAnderson order parameter. In the RSB picture, P J ( q ) consists of a countable infinity of δ functions of varying weights densely filling the range between ± q EA while P q ( ) is a smooth function between ± q EA with delta functions at ± q EA . Thus, one can, in principle, distinguish between the two classes of pictures by examining P q ( ) near the origin (and thus away from q EA ). A measure of the weight near the origin of the overlap distribution is I J ( q ),

$$I _ { \mathcal { J } } ( q ) = \int _ { - q _ { 0 } } ^ { + q _ { 0 } } d q P _ { \mathcal { J } } ( q ). \quad \quad ( 9 ) \quad \text{TA} \quad \text{diff}$$

We refer to the disorder average of I J ( q ) as I ( q ). The choice q 0 = 0 2 . has been used in many past studies [38, 44] to distinguish the RSB and two-state pictures and, in this work, we investigate the statistics of I J ( q 0 = 0 2). . In the following, we use the symbols I J and I L as abbreviations for I J (0 2) and its disorder av-. erage I (0 2) for size . L , respectively.

## V. RESULTS

In this section, we present results for the spin stiffness (VA), the order parameter distribution near zero, I J (VB) and the correlation of I J and λ J (VC). The main result of this section is that λ J increases with system size and that stiff samples have small values of I J .

## A. Spin stiffness

Figure 1 shows the free-energy difference or, for T = 0, ground-state energy difference ∆ F between TBC and PBC for temperatures T = 0, 0 2, . and 0 42 as a func-. tion of system size L . The straight lines are best fits to the functional form ∆ F ∼ aL θ . The fits for θ are shown in Table II. The result θ T ( = 0) = 0 197 . ± 0 017 . is in reasonable agreement though at the low end of previous measurements of θ carried out at zero temperature [37, 66-68]. Note that the stiffness exponent has not previously been measured at nonzero temperature. We see that θ decreases as temperature increases. Presumably, this is a finite-size effect because θ is expected to have a single asymptotic value throughout the low-temperature phase [30].

Next consider the sample stiffness measure λ J , defined in Eq. (7). Let G L ( λ ) be the cumulative distribution

FIG. 1: Free-energy change ∆ F vs system size L for T = 0, 0 2, . and 0 42. . The straight lines are fits of the form ∆ F ∼ aL θ .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0438v2-figure-9709d826-da92-4953-bd0a-4ab0932b42f2.jpg)

TABLE II: Estimates of the stiffness exponents θ and θ λ for different temperatures T .

| T   | 0           | 0 . 2       | 0 . 42      |
|-----|-------------|-------------|-------------|
| θ   | 0 . 197(17) | 0 . 189(17) | 0 . 169(12) |
| θ λ | -           | 0 . 290(30) | 0 . 268(20) |

function for λ . The left panel of Fig. 2 is a log plot of 1 -G L ( λ ), the complementary cumulative distribution function, for λ at T = 0 42, and sizes 4 through 12. . The nearly straight line behavior of log(1 -G L ( λ )) is indicative of a nearly exponential tail and suggests a data collapse if λ is scaled by a characteristic λ char ( L ) given by the slope of the line. Since the tail is not perfectly straight, we instead define λ char ( L ) in terms of medianlike quantities. If the distribution were exactly exponential then 1 -G λ ( ) = e -λ/λ char and 1 -G λ ( char log b ) = 1 /b for any b . If the distribution is not perfectly exponential, λ char depends on b so there is some ambiguity in the definition. We choose b such that λ char is obtained from the tail of the distribution but not so far into the tail that the statistics are poor. For T = 0 2 we choose . b = 2 so λ char is defined as the median divided by log 2. For T = 0 42 we choose . b = 10 to ensure that λ char is obtained from the tail of the distribution. The right panels of Figs. 2 and 3 show 1 -G L ( λ/λ char ( L )) for T = 0 42 . and T = 0 2, respectively, and reveal that all of the cu-. mulative distributions collapse onto the same curve when scaled by λ char ( L ).

Figure 4 shows λ char ( L ) vs log L for T = 0 2 and . T = 0 42. . Since λ char ( L ) is a stiffness measure, we can extract a new stiffness exponent θ λ from a fit to the form,

$$\lambda _ { \text{char} } ( L ) \sim a L ^ { \theta _ { \lambda } }. \quad \quad \quad ( 1 0 )$$

The values of θ λ , given in Table II, are larger than θ obtained from the average free energy difference but close to the value, 0 27, found in Ref. [67] using aspect ratio .

scaling. Presumably, the asymptotic values of θ and θ λ are the same. We prefer the larger value, θ λ because it is obtained from the tail of the stiffness distribution so we believe it reflects the large-size behavior more accurately than the average free energy difference that defines θ . Aspect ratio scaling is an independent way to minimize finite-size effects and it is interesting that these two approaches yield the same answer within error bars.

It seems clear that λ char ( L ) →∞ as L →∞ . At least on a coarse scale, the full distribution G L ( λ ) also scales with λ char ( L ). A closer look at G near the head of the distribution for T = 0 42 shows that the data collapse is . not perfect and there are significant finite-size corrections near λ = 0. Figure 5 shows 1 -G L ( λ/λ char ( L )) vs λ in the region near λ = 0 for T = 0 42. . Note that curves do not collapse perfectly and that 1 -G L ( λ/λ char ( L )) appears to be increasing with L . Figure 6 is the same plot for T = 0 2. Because . G L (0) is so small for T = 0 2, the error . bars are too large to discern whether there is a trend with L . A reasonable hypothesis is that there is an asymptotic L →∞ scaling function G ∞ ( z ) where z = λ/λ char such that G L ( λ ) → G ∞ ( λ/λ char ( L )). The straight line behavior of log(1 -G ∞ ( z )) for large z and increasing trend with L for small z suggests that G ∞ (0) = 0 and G ∞ ( z ) is exponential for z glyph[greatermuch] 1. In more physical terms, if G ∞ ( z ) exists and is zero for z → 0 + , it means that a single boundary condition dominates the TBC ensemble almost surely, i.e., the dominant boundary condition almost always has a much lower free energy than the other seven boundary conditions. A more complicated possibility is that G ∞ (0 + ) &gt; 0. The consequences of these possibilities for the RSB vs two-state pictures are discussed in Sec. VI.

It is noteworthy that for the system sizes accessible to simulations, λ char ( L ) is sufficiently small that the TBC ensemble contains a mixture of several competing boundary conditions for a substantial fraction of samples. The disorder average I L is dominated by these samples and is therefore not characteristic of the largeL behavior when λ char ( L ) is expected to be large. In what follows we circumvent this difficulty by extrapolating first in λ and then in L , making use of the fact our data contain a relatively large dynamic range in λ .

## B. Order parameter near q = 0

Figure 7 shows I L , the disorder average of the integrated order parameter distribution with | q | &lt; 0 2, as a . function of size L for temperatures T = 0 2 and 0 42, as . . well as for both PBC and TBC. For both boundary conditions, I L is, within error bars, independent of system size. The results for PBC are identical to those obtained using parallel tempering Monte Carlo [38, 44]. In fact, in Fig. 8 we show a scatter plot of I J computed both with population annealing and parallel tempering Monte Carlo. The data are strongly correlated and both methods yield the same results within statistical errors.

FIG. 2: (Left panel) Linear-log plot of 1 -G L ( λ ) (the complementary cumulative distribution function) vs λ for sizes L = 4 through 12 at T = 0 42. (Right panel) 1 . -G L ( λ/λ char ( L )) vs λ/λ char ( L ).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0438v2-figure-e4d885de-7fd6-44c2-9a38-3c3315951f52.jpg)

FIG. 3: (Left panel) Linear-log plot of 1 -G L ( λ ) (the complementary cumulative distribution function) vs λ for sizes L = 4 through 10 at T = 0 2. . (Right panel) 1 -G L ( λ/λ char ( L )) vs λ/λ char ( L ).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0438v2-figure-193538f2-e807-4cc8-9609-47388acb7cec.jpg)

The constancy of I L has been taken as strong evidence for the RSB picture because the two-state picture predicts I L should decrease as L -θ . However, in what follows we argue that in TBC ultimately I L → 0 for very large L . On first glance the results for TBC are surprising since I TBC L is larger by more than a factor of two than I PBC L . The explanation is that for many samples the TBC ensemble contains several boundary conditions with significant weight and the overlap between spin configurations with different boundary conditions will tend to have small values of q due to the existence of a relative domain wall. We shall return to this important point in Sec. VI.

## C. Order parameter near q = 0 vs spin stiffness

Figure 9 is a scatter plot showing many of the disorder samples at T = 0 42 (left panel) and . T = 0 2 (right panel) . for all the sizes studied using TBC. Each point on the plot represents a sample J . The x coordinate of the point is λ J = log [ f J / (1 -f J )] and the y coordinate is I J .

FIG. 4: Log-log plot of λ char ( L ) vs L for T = 0 2 and . T = 0 42. . The straight lines represent fits of the form λ char ( L ) ∼ aL θ λ .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0438v2-figure-c967a276-8795-4ff0-a5c1-b2339ab5cbc4.jpg)

FIG. 5: (Left panel) 1 -G L ( λ ) vs λ for system sizes L = 4 through 12 at T = 0 42 in the region near . λ = 0. (Right panel) 1 -G L ( λ/λ char ( L )) vs λ/λ char ( L ). Note that 1 -G L (0) increases slowly with L .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0438v2-figure-3605da24-775c-4c7a-af04-b619948883e5.jpg)

Figure 10 is the same as Fig. 9 but with each system size on a separate plot for T = 0 2. The qualitative features of . the plots are the same for each size although, as described above, the distribution of λ 's shifts to larger values for larger L . These figures together with the behavior of the λ distribution constitute the main results of our paper and motivate our conclusion that I J → 0 almost surely as L →∞ in thermal boundary conditions.

Samples J for which I J or 1 -f J are exactly zero within the precision of the simulations are not shown on these log-log plots since log I J = -∞ or λ J = ∞ . The fractions of such samples are given in Table III. Note that actual values of I J or 1 -f J are are never exactly zero for finite systems; zeros correspond to values smaller than can be represented by the finite population sizes used in the simulations. It is important to note that the trends shown in Fig. 9 continue to hold for the large values of λ that are omitted from this figure. Including all sizes, there are 216 samples with 1 -f J = 0 for T = 0 42 . and 2996 such samples for T = 0 2. . Of these, only seven

FIG. 6: (Left panel) 1 -G L ( λ ) vs λ for system sizes L = 4 through 10 at T = 0 2 in the region near . λ = 0. (Right panel) 1 -G L ( λ/λ char ( L )) vs λ/λ char ( L ).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0438v2-figure-ea9428a6-1733-4a2e-8462-1c600d31030a.jpg)

FIG. 7: I L vs L for PBC and TBC at temperature T = 0 42 (left . panel) and T = 0 2 (right panel). . The data seem independent of system size, suggesting an RSB interpretation of the data.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0438v2-figure-6aaed808-fbc5-4fd1-8505-cbfaa2b3586a.jpg)

samples for T = 0 42 and 38 for . T = 0 2 are measured . to have nonzero values of I J . The average value of I for only those samples with 1 -f J = 0 are 2 5 . × 10 -8 and 10 -4 for T = 0 42 and . T = 0 2, respectively. .

A striking feature of Figs. 9 and 10 is that there is a bounding curve that becomes a straight line for large λ with most samples lying below that curve. Why are there two classes of samples, with most samples below the curve and a few above it? We speculate that the samples below the curve have nonzero values of I J as a result of the overlap between spin configurations with different boundary conditions. The contribution to the overlap between different boundary conditions in the TBC ensemble cannot exceed 2 f (1 -f ) in the limit f → 1, so that if this is the primary mechanism producing small overlaps in sample J then log I J &lt; ( -λ J + log 2). The straight lines in Fig. 9 are defined by log I = ( -λ +log 2). Thus, for most samples, we believe that the primary contribution to I J comes from the overlap between different boundary conditions. For the rare samples above the bounding curve,

FIG. 8: Log-log scatter plot of I J for PBC for each sample at temperature T = 0 20 and . L = 10 computed with parallel tempering Monte Carlo (PT, vertical axis) and population annealing Monte Carlo (PA, horizontal line). The data show both methods yields the same result for each sample within statistical errors.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0438v2-figure-052e69be-15a6-409b-ba87-7bd2a5affb21.jpg)

TABLE III: Fraction of samples with I J = 0 and f J = 1 for different sizes, temperatures, and boundary conditions.

| PBC                               | PBC   | PBC   | PBC   | PBC   | PBC   |
|-----------------------------------|-------|-------|-------|-------|-------|
| System size L                     | 4     | 6     | 8     | 10    | 12    |
| Fraction I J = 0 ( T = 0 . 42)    | 0.21  | 0.19  | 0.16  | 0.16  | 0.19  |
| Fraction of I J = 0 ( T = 0 . 2)  | 0.60  | 0.57  | 0.55  | 0.54  | -     |
| Fraction of f J = 1 ( T = 0 . 42) | -     | -     | -     | -     | -     |
| Fraction of f J = 1 ( T = 0 . 2)  | -     | -     | -     | -     | -     |

| TBC                               |   TBC |   TBC |   TBC |    TBC | TBC   |
|-----------------------------------|-------|-------|-------|--------|-------|
| System size L                     | 4     | 6     | 8     | 10     | 12    |
| Fraction I J = 0 ( T = 0 . 42)    | 0.05  | 0.04  | 0.03  |  0.03  | 0.03  |
| Fraction of I J = 0 ( T = 0 . 2)  | 0.35  | 0.33  | 0.3   |  0.28  | -     |
| Fraction of f J = 1 ( T = 0 . 42) | 0.009 | 0.009 | 0.008 |  0.008 | 0.010 |
| Fraction of f J = 1 ( T = 0 . 2)  | 0.15  | 0.16  | 0.15  |  0.15  | -     |

the primary contribution to I J must come from small overlaps within the dominant boundary condition.

Asecond important feature of Figs. 9 and 10 is that the rare samples above the bounding curve have I J roughly uniformly distributed on a logarithmic scale between the bounding curve and 1. On a linear scale this means that for large λ almost all of these samples have small values of I J . Let ρ x y ( | ) be the conditional probability density for x = I J conditioned on y = λ J . If this distribution is exactly log uniform above the bounding line then the part of the distribution above the line would take the

FIG. 9: (Color online) Scatter plots showing all disorder realizations for all system sizes at T = 0 42 (left . panel) and T = 0 2 (right . panel). Each point represents a sample J located at x coordinate λ J and y coordinate I J . Red diamonds represent L = 4, blue crosses L = 6, green squares L = 8, purple triangles L = 10, and orange plus symbols L = 12.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0438v2-figure-ae251e4e-daa4-4bda-b5af-47844672c688.jpg)

FIG. 10: (Color online) Same as Fig. 9 but for each system size in a separate panel and T = 0 2. . Again, red diamonds represent L = 4, blue crosses L = 6, green squares L = 8, and purple triangles L = 10.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0438v2-figure-51240ea5-c912-4c9e-a4bd-ff149019c995.jpg)

form,

$$- \quad \rho ( x | y ) = \frac { 1 - w ( y ) } { x y } \quad \text{for} \quad x > 2 \exp ( - y ), \ \ ( 1 1 )$$

where w y ( ) is the fraction samples below the bounding line. Figure 11 shows histograms of I J values of the samples of all sizes that lie above the bounding line for the two temperatures. The position α along the horizontal axis is the scaled logarithmic distance between the bounding line and one. That is, α J = -log I J / λ [ -log(2)] so that zero corresponds to large values, I J ≈ 1 while one corresponds to I J on the bounding line. For T = 0 2 the distribution is indeed relatively uniform on . a logarithmic scale while for T = 0 42 it is skewed to a . small value of I .

)

α

(

N

30

20

10

0

0

0.2

0.4

0.6

0.8

1

α

FIG. 11: Histogram N ( α ) for α J = log( I J ) / [ -λ J +log(2)] for T = 0 42 (left panel) and . T = 0 2 (right panel). . α = 1 corresponds to the small values of I at the bounding line.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0438v2-figure-5170dc34-f245-41ae-8bb6-65674d6e2f51.jpg)

## VI. DISCUSSION

Three salient features of the data are apparent from Figs. 4, 9, and 11.

- I. Typical values of the sample stiffness λ J increase with system sizes L , as described by λ char ( L ).
- II. Most samples have I J less than a bounding curve described by 2 e -λ for large λ .
- III. Samples with I J above the bounding curve have I J distributed more or less uniformly in log I between the bounding curve and one.

We conjecture that these features hold for arbitrarily large L and all temperatures in the low-temperature phase. Assuming the above statements are asymptotically correct, we can draw some strong conclusions about how I L behaves for sufficiently large L that λ char ( L ) glyph[greatermuch] 1. Note that these system sizes are far larger than are accessible in our simulations but, given the large dynamic range in λ we can extrapolate to these sizes by first extrapolating in λ . For large L λ , J is nearly always large according to Statement (I). Furthermore, due to Statements (II) and (III), I J is, almost always small when λ J is large. Thus I J is nearly always small when L is large. This conclusion is the main result of our analysis. It is consistent with two-state pictures but not consistent with the RSB picture.

In addition to being consistent with our data, these conjectures are quite plausible. Statement (I) asserts that λ J is a measure of sample stiffness and that in the low-temperature phase, almost all samples become stiff for large system sizes. As discussed above, Statement (II) asserts that in TBC large values of I J arise mostly from the overlap between two different boundary conditions. Statement (III) asserts that the free-energy cost of a large excitation in the dominant boundary condition is more or less uniformly distributed between λ/β and 0.

We can make the arguments more quantitative using a simple model of how the disorder average I L will behave for large L . Let ρ x y ( | ) be the conditional probability density for I J = x conditioned on λ J = y . Based on statements (II) and (III) we propose the form,

$$\mathfrak { P } \, \Big | \quad & \rho ( x | y ) = w ( y ) \delta [ x, 2 \exp ( - y ) ] + \frac { 1 - w ( y ) } { x y } \theta [ x - 2 \exp ( - y ) ], \\ & \quad \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \ dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots.$$

where w y ( ) is the fraction of samples at fixed λ below the bounding curve, and θ x ( ) and δ x, y ( ) are the Heaviside function and the δ -function, respectively. The first term conservatively places all the samples below the bounding curve on the curve itself. The second term represents the samples above the bounding curve with the conservative assumption that the distribution of I J above the curve is uniform on a log scale as in Eq. (11). For purposes of the following calculation we assume that w is a constant independent of y but the qualitative conclusions do not depend on this assumption. Finally, we assume that the distribution of λ obeys a size independent form G ∞ ( z ) for the scaled variable z = λ/λ char ( L ). Note that we have assumed that the dependence of I L on L is entirely through λ char ( L ) and that the conditional probability ρ x y ( | ) is independent of L . These assumptions yield an explicit expression for I L as a function of λ char ( L ),

$$\ r v e \quad & I _ { L } = \frac { 1 } { \lambda _ { \text{char} } ( L ) } \int _ { 0 } ^ { \infty } d G _ { \infty } ( z ) \int _ { 0 } ^ { \infty } d x \ x \rho ( x | z \lambda _ { \text{char} } ( L ) ). \\ I _ { \tau } \quad & \text{d} \dots \dots \colon \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon = \colon$$

Plugging in the ansatz of Eq. (12) and an exponential form for the scaled λ distribution, 1 -G ∞ ( z ) = e -z , yields a somewhat complex expression involving exponential integrals whose asymptotic large λ behavior is,

$$\overset { \circ } { \cdot } \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd. \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd. \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cdk \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cdj \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cdjs \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cdJ \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd3 \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd1 \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \cd� \\ \\
 \\$$

Using Eq. (10) and assuming that asymptotically θ λ = , θ we recover the prediction of the two-state picture that I L ∼ L -θ , however, with a logarithmic correction that arises from the assumption of a log-uniform distribution for I J .

While the above assumptions lead to an explicit asymptotic expression for I L as a function L , this expression should not be taken too seriously. However, the qualitative conclusion that I J → 0 for almost all samples as L →∞ is robust and depends only on the asymptotic validity of statements (I) - (III) above.

Given that λ char ( L ) is increasing with L and that I decreases with increasing λ , why is I L nearly constant for the sizes studied in our TBC simulations of the 3D EA model? We believe this conundrum can be explained, at least in part, by the fact that the main contribution to I L comes from samples with small values of λ . For T = 0 2, . more than half the contribution to I L comes from samples with λ &lt; 1 and more than 80% from λ &lt; 2,

and these fractions are even higher for T = 0 42. . The head of the λ distribution, has very little dependence on L , as can be seen in Figs. 2 and 3. Furthermore, the bounding curves in Fig. 9 are nearly flat in the small λ region. Thus several effects come into play in keeping I L nearly independent of L . First, the main contribution to I L is from samples with small stiffness. Second, the fraction of samples with small stiffness decreases by only a small amount for the sizes studied and, finally, I does not depend much on λ for small λ . One would have to go to much larger sizes before I L would decrease according to the predicted asymptotic power law L -θ .

In the foregoing, we have assumed that G L (0) → 0 as L → ∞ or, equivalently, if G ∞ ( z ) exists, G ∞ (0 + ) = 0. We now consider the consequences of an alternate assumption that λ char ( L ) → ∞ and G ∞ ( z ) exists but G ∞ (0 + ) &gt; 0. This possibility cannot be ruled out by the data although if it holds, it appears that G ∞ (0 + ) is quite small. In physical terms G ∞ (0 + ) &gt; 0 means that even for very large sizes, a fraction G ∞ (0 + ) of samples has a mixed ensemble of boundary conditions in TBC while the remaining samples have only a single boundary condition in the TBC ensemble. This scenario would imply that the 3D EA model in TBC is divided into two classes of disorder realizations, one of which, with weight (1 -G ∞ (0 + )), has I L = 0 and the other, with weight G ∞ (0 + ), has I L &gt; 0. This possibility seems unlikely but is not contradicted by the data. It has no straightforward explanation in either two-state or RSB pictures.

Our hypothesis is that thermal boundary conditions and periodic boundary conditions have the same behavior in the limit of large system sizes. We use thermal boundary conditions as a tool to improve the extrapolation to large system sizes from the very small system sizes accessible in simulations. It is known that coupling dependent boundary conditions are not equivalent to periodic boundary conditions and are not suitable for understanding properties of the spin glass phase because they could be used to select a single pure state even if coupling independent boundary conditions admit many pure states. The status of thermal boundary conditions with regard to coupling dependence is not clear. On the one hand, the TBC ensemble contains different mixtures of boundary conditions for different choices of couplings. On the other hand, the particular mixture of the 8 boundary conditions is chosen by the system itself and is not imposed externally. As discussed in Sec. II B, our intuition is that TBC minimizes finite size effects rather than introducing spurious physics but this question requires further investigation. In any case, we have provided compelling evidence that the 3D EA model in thermal boundary conditions is best described by a picture with a single pair of pure states in each finite volume.

## VII. SUMMARY

We have introduced two techniques with the aim of extrapolating to the large system-size behavior of finitedimensional spin glasses at low temperature. First, we use thermal boundary conditions to minimize the effect of domain walls induced by boundary conditions. Second, we use a natural measure of sample stiffness defined within thermal boundary conditions and extrapolate to large values of the sample stiffness. By noting that the sample stiffness increases with system size we then obtain an extrapolation in system size. The dynamic range in sample stiffness in the data is sufficiently large that a qualitative extrapolation is readily apparent. The conclusion is that nearly all large samples will have essentially no weight in the overlap distribution near zero overlap. The analysis also explains why this qualitative behavior cannot be seen using a direct extrapolation in system size for the small sizes studied. Our conclusions are consistent with two-state pictures but are inconsistent with the mean field, replica symmetry breaking picture.

Our results hold for thermal boundary conditions. We believe that thermal boundary conditions are equivalent to other coupling independent boundary conditions so that our conclusions about the infinite volume limit also apply to the more familiar periodic boundary conditions. However, it is important to investigate the equivalence of thermal and periodic boundary conditions.

Our numerical simulations used population annealing. This Monte Carlo algorithm has not been used before for large-scale studies in statistical physics. We found that it is an effective computational tool with several advantages over parallel tempering, the standard computational method in the field. We believe that population annealing will be useful for other hard problems in statistical physics and related fields. Similarly, thermal boundary conditions and extrapolating in sample stiffness are general methods that should be useful in studying other finite-dimensional disordered systems.

## Acknowledgments

H.G.K. acknowledges support from the NSF (Grant No. DMR-1151387) and would like to thank Georg Schneider &amp; Sohn for creating TAP4. J.M. and W.W. acknowledge support from NSF (Grant No. DMR1208046). We are grateful to Alan Middleton, Dan Stein, and Chuck Newman for reading an early version of the manuscript and for useful discussions. We acknowledge the contribution of Burcu Yucesoy in providing comparison data from parallel tempering simulations and for useful discussions. We thank the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources (Stampede Cluster), ETH Zurich for CPU time on the Brutus and Euler clusters, and Texas A&amp;M University for access to their Eos and Lonestar clusters. We especially thank O. Byrde for beta

- [1] K. Binder and A. P. Young, Spin glasses: Experimental facts, theoretical concepts and open questions , Rev. Mod. Phys. 58 , 801 (1986).
- [2] M. M´ ezard, G. Parisi, and M. A. Virasoro, Spin Glass Theory and Beyond (World Scientific, Singapore, 1987).
- [3] A. P. Young, ed., Spin Glasses and Random Fields (World Scientific, Singapore, 1998).
- [4] A. K. Hartmann and H. Rieger, New Optimization Algorithms in Physics (Wiley-VCH, Berlin, 2004).
- [5] D. L. Stein and C. M. Newman, Spin Glasses and Complexity , Primers in Complex Systems (Princeton University Press, 2013).
- [6] A. Lucas, Ising formulations of many NP problems , Front. Physics 12 , 5 (2014).
- [7] M. H. Freedman, A. Y. Kitaev, M. J. Larsen, and Z. Wang, Topological Quantum Computation (2002), (arXiv:quant-ph/0101025).
- [8] A. Y. Kitaev, Fault-tolerant quantum computation by anyons , Ann. Phys. 303 , 2 (2003).
- [9] H. Bombin and M. A. Martin-Delgado, Topological Quantum Distillation , Phys. Rev. Lett. 97 , 180501 (2006).
- [10] D. Nigg, M. Mueller, E. A. Martinez, P. Schindler, M. Hennrich, T. Monz, M. A. Martin-Delgado, and R. Blatt, Experimental Quantum Computations on a Topologically Encoded Qubit , Science 345 , 302 (2014).
- [11] E. Dennis, A. Kitaev, A. Landahl, and J. Preskill, Topological quantum memory , J. Math. Phys. 43 , 4452 (2002).
- [12] H. G. Katzgraber, H. Bombin, and M. A. MartinDelgado, Error Threshold for Color Codes and Random 3-Body Ising Models , Phys. Rev. Lett. 103 , 090501 (2009).
- [13] H. Bombin, R. S. Andrist, M. Ohzeki, H. G. Katzgraber, and M. A. Martin-Delgado, Strong Resilience of Topological Codes to Depolarization , Phys. Rev. X 2 , 021004 (2012).
- [14] N. G. Dickson, M. W. Johnson, M. H. Amin, R. Harris, F. Altomare, A. J. Berkley, P. Bunyk, J. Cai, E. M. Chapple, P. Chavez, et al., Thermally assisted quantum annealing of a 16-qubit problem , Nat. Comm. 4 , 1903 (2013).
- [15] S. Boixo, T. Albash, F. M. Spedalieri, N. Chancellor, and D. A. Lidar, Experimental signature of programmable quantum annealing , Nat. Comm. 4 , 2067 (2013).
- [16] H. G. Katzgraber, F. Hamze, and R. S. Andrist, Glassy Chimeras Could Be Blind to Quantum Speedup: Designing Better Benchmarks for Quantum Annealing Machines , Phys. Rev. X 4 , 021008 (2014).
- [17] T. F. Rønnow, Z. Wang, J. Job, S. Boixo, S. V. Isakov, D. Wecker, J. M. Martinis, D. A. Lidar, and M. Troyer, Defining and detecting quantum speedup , Science 345 (2014).
- [18] G. Parisi, Infinite number of order parameters for spinglasses , Phys. Rev. Lett. 43 , 1754 (1979).
- [19] G. Parisi, The order parameter for spin glasses: a function on the interval 0 -1, J. Phys. A 13 , 1101 (1980).
- [20] G. Parisi, Order parameter for spin-glasses , Phys. Rev. Lett. 50 , 1946 (1983).
- [21] R. Rammal, G. Toulouse, and M. A. Virasoro, Ultra-
- metricity for physicists , Rev. Mod. Phys. 58 , 765 (1986).
- [22] S. F. Edwards and P. W. Anderson, Theory of spin glasses , J. Phys. F: Met. Phys. 5 , 965 (1975).
- [23] G. Parisi, Some considerations of finite dimensional spin glasses , J. Phys. A 41 , 324002 (2008).
- [24] D. Sherrington and S. Kirkpatrick, Solvable model of a spin glass , Phys. Rev. Lett. 35 , 1792 (1975).
- [25] D. Panchenko, The Sherrington-Kirkpatrick model: An overview , J. Stat. Phys. 149 , 362 (2012).
- [26] W. L. McMillan, Domain-wall renormalization-group study of the two-dimensional random Ising model , Phys. Rev. B 29 , 4026 (1984).
- [27] A. J. Bray and M. A. Moore, Scaling theory of the ordered phase of spin glasses , in Heidelberg Colloquium on Glassy Dynamics and Optimization , edited by L. Van Hemmen and I. Morgenstern (Springer, New York, 1986), p. 121.
- [28] D. S. Fisher and D. A. Huse, Ordered phase of short-range Ising spin-glasses , Phys. Rev. Lett. 56 , 1601 (1986).
- [29] D. S. Fisher and D. A. Huse, Absence of many states in realistic spin glasses , J. Phys. A 20 , L1005 (1987).
- [30] D. S. Fisher and D. A. Huse, Equilibrium behavior of the spin-glass ordered phase , Phys. Rev. B 38 , 386 (1988).
- [31] C. M. Newman and D. L. Stein, Multiple states and thermodynamic limits in short-ranged Ising spin-glass models , Phys. Rev. B 46 , 973 (1992).
- [32] C. M. Newman and D. L. Stein, Non-mean-field behavior of realistic spin glasses , Phys. Rev. Lett. 76 , 515 (1996).
- [33] C. M. Newman and D. L. Stein, Simplicity of state and overlap structure in finite-volume realistic spin glasses , Phys. Rev. E 57 , 1356 (1998).
- [34] C. M. Newman and D. L. Stein, The state(s) of replica symmetry breaking: Mean field theories vs short-ranged spin glasses , J. Stat. Phys. 106 , 213 (2002).
- [35] C. M. Newman and D. L. Stein, TOPICAL REVIEW: Ordering and broken symmetry in short-ranged spin glasses , J. Phys.: Condensed Matter 15 , 1319 (2003).
- [36] N. Read, Short-range Ising spin glasses: the metastate interpretation of replica symmetry breaking , Phys. Rev. E 90 , 032142 (2014).
- [37] M. Palassini and A. P. Young, Triviality of the ground state structure in Ising spin glasses , Phys. Rev. Lett. 83 , 5126 (1999).
- [38] H. G. Katzgraber, M. Palassini, and A. P. Young, Monte Carlo simulations of spin glasses at low temperatures , Phys. Rev. B 63 , 184422 (2001).
- [39] H. G. Katzgraber and A. P. Young, Monte Carlo simulations of spin-glasses at low temperatures: Effects of free boundary conditions , Phys. Rev. B 65 , 214402 (2002).
- [40] R. Alvarez Ba˜os, A. Cruz, L. A. Fernandez, J. M. Giln Narvion, A. Gordillo-Guerrero, M. Guidetti, A. Maiorano, F. Mantovani, E. Marinari, V. Martin-Mayor, et al., Nature of the spin-glass phase at experimental length scales , J. Stat. Mech. P06026 (2010).
- [41] J. Houdayer, F. Krzakala, and O. C. Martin, Large-scale low-energy excitations in 3-d spin glasses , Eur. Phys. J. B. 18 , 467 (2000).
- [42] F. Krzakala and O. C. Martin, Spin and link overlaps in 3-dimensional spin glasses , Phys. Rev. Lett. 85 , 3013

(2000).

- [43] H. G. Katzgraber and A. P. Young, Monte Carlo studies of the one-dimensional Ising spin glass with power-law interactions , Phys. Rev. B 67 , 134410 (2003).
- [44] B. Yucesoy, H. G. Katzgraber, and J. Machta, Evidence of Non-Mean-Field-Like Low-Temperature Behavior in the Edwards-Anderson Spin-Glass Model , Phys. Rev. Lett. 109 , 177204 (2012).
- [45] A. Billoire, A. Maiorano, E. Marinari, V. Martin-Mayor, and D. Yllanes, The cumulative overlap distribution function in realistic spin glasses , Phys. Rev. B 90 , 094201 (2014).
- [46] M. Palassini and A. P. Young, Nature of the spin glass state , Phys. Rev. Lett. 85 , 3017 (2000).
- [47] A. Billoire, L. A. Fernandez, A. Maiorano, E. Marinari, V. Martin-Mayor, G. Parisi, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, and D. Yllanes, Comment on 'Evidence of Non-Mean-Field-Like Low-Temperature Behavior in the Edwards-Anderson Spin-Glass Model' , Phys. Rev. Lett. 110 , 219701 (2013).
- [48] B. Yucesoy, H. G. Katzgraber, and J. Machta, Yucesoy, Katzgraber, and Machta reply: , Phys. Rev. Lett. 110 , 219702 (2013).
- [49] A. A. Middleton, Extracting thermodynamic behavior of spin glasses from the overlap function , Phys. Rev. B 87 , 220201 (2013).
- [50] C. Monthus and T. Garel, Typical versus averaged overlap distribution in spin glasses: Evidence for droplet scaling theory , Phys. Rev. B 88 , 134204 (2013).
- [51] A. Billoire, A. Maiorano, E. Marinari, V. Martin-Mayor, and D. Yllanes, The cumulative overlap distribution function in realistic spin glasses , ArXiv e-prints (2014), 1406.1639.
- [52] J. W. Landry and S. N. Coppersmith, Ground states of two-dimensional ± J Edwards-Anderson spin glasses , Phys. Rev. B 65 , 134404 (2002).
- [53] C. K. Thomas and A. A. Middleton, Matching Kasteleyn cities for spin glass ground states , Phys. Rev. B 76 , 220406(R) (2007).
- [54] K. Hukushima, Domain-wall free energy of spin-glass models: Numerical method and boundary conditions , Phys. Rev. E 60 , 3606 (1999).
- [55] M. Sasaki, K. Hukushima, H. Yoshino, and H. Takayama, Temperature Chaos and Bond Chaos in EdwardsAnderson Ising Spin Glasses: Domain-Wall Free-Energy
- Measurements , Phys. Rev. Lett. 95 , 267203 (2005).
- [56] M. Sasaki, K. Hukushima, H. Yoshino, and H. Takayama, Scaling Analysis of Domain-Wall Free Energy in the Edwards-Anderson Ising Spin Glass in a Magnetic Field , Phys. Rev. Lett. 99 , 137202 (2007).
- [57] M. Hasenbusch, Monte-Carlo simulation with fluctuating boundary-conditions , Physica A 197 , 423 (1993).
- [58] K. Hukushima and Y. Iba, in THE MONTE CARLO METHOD IN THE PHYSICAL SCIENCES: Celebrating the 50th Anniversary of the Metropolis Algorithm , edited by J. E. Gubernatis (AIP, 2003), vol. 690, pp. 200-206.
- [59] J. Machta, Population annealing with weighted averages: A Monte Carlo method for rough free-energy landscapes , Phys. Rev. E 82 , 026704 (2010).
- [60] J. Machta and R. Ellis, Monte Carlo methods for rough free energy landscapes: Population annealing and parallel tempering , J. Stat. Phys. 144 , 541 (2011).
- [61] K. Hukushima and K. Nemoto, Exchange Monte Carlo method and application to spin glass simulations , J. Phys. Soc. Jpn. 65 , 1604 (1996).
- [62] A. Doucet, N. de Freitas, and N. Gordon, eds., Sequential Monte Carlo Methods in Practice (Springer-Verlag, New York, 2001).
- [63] B. Yucesoy, J. Machta, and H. G. Katzgraber, Correlations between the dynamics of parallel tempering and the free-energy landscape in spin glasses , Phys. Rev. E 87 , 012104 (2013).
- [64] W. Wang, J. Machta, and H. G. Katzgraber, Population annealing for large scale spin glass simulations (2014), in preparation.
- [65] W. Wang, J. Machta, and H. G. Katzgraber, Finding ground states of Ising spin glasses using population annealing Monte Carlo (2014), in preparation.
- [66] A. K. Hartmann, Scaling of stiffness energy for threedimensional ± J Ising spin glasses , Phys. Rev. E 59 , 84 (1999).
- [67] A. C. Carter, A. J. Bray, and M. A. Moore, Aspect-ratio scaling and the stiffness exponent θ for Ising spin glasses , Phys. Rev. Lett. 88 , 077201 (2002).
- [68] S. Boettcher, Stiffness exponents for lattice spin glasses in dimensions d = 3 , . . . , 6, The European Physical Journal B - Condensed Matter and Complex Systems 38 , 83 (2004).
Wenlong Wang Jonathan Machta Helmut G. Katzgraber
Evidence against a mean field description of short-range spin glasses revealed through thermal boundary conditions
A theoretical description of the low-temperature phase of short-range spin
glasses has remained elusive for decades. In particular, it is unclear if
theories that assert a single pair of pure states, or theories that are based
infinitely many pure states-such as replica symmetry breaking-best describe
realistic short-range systems. To resolve this controversy, the
three-dimensional Edwards-Anderson Ising spin glass in thermal boundary
conditions is studied numerically using population annealing Monte Carlo. In
thermal boundary conditions all eight combinations of periodic vs antiperiodic
boundary conditions in the three spatial directions appear in the ensemble with
their respective Boltzmann weights, thus minimizing finite-size corrections due
to domain walls. From the relative weighting of the eight boundary conditions
for each disorder instance a sample stiffness is defined, and its typical value
is shown to grow with system size according to a stiffness exponent. An
extrapolation to the large-system-size limit is in agreement with a description
that supports the droplet picture and other theories that assert a single pair
of pure states. The results are, however, incompatible with the mean-field
replica symmetry breaking picture, thus highlighting the need to go beyond
mean-field descriptions to accurately describe short-range spin-glass systems.