{'fraction_non_alphanumeric': 0.03097201714265135, 'fraction_numerical': 0.012929016458385853, 'mean_word_length': 4.528170416086004, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 7, '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}}
Sampling the ground-state magnetization of d-dimensional p-body Ising models

We demonstrate that a recently introduced heuristic optimization algorithm [Phys. Rev. E 83, 046709 (2011)] that combines a local search with triadic crossover genetic updates is capable of sampling nearly uniformly among ground-state configurations in spin-glass-like Hamiltonians with p-spin interactions in d space dimensions that have highly degenerate ground states. Using this algorithm we probe the zero-temperature ferromagnet to spin-glass transition point q_c of two example models, the disordered version of the two-dimensional three-spin Baxter-Wu model [q_c = 0.1072(1)] and the three-dimensional Edwards-Anderson model [q_c = 0.2253(7)], by computing the Binder ratio of the ground-state magnetization.


I. INTRODUCTION
Disordered systems with degenerate ground states often have an exponentially-large 1 number of states that minimize the Hamiltonian (cost function) of the problem. A hallmark model system to study degenerate ground states is the Edwards Anderson Ising spin-glass model with bimodal (±J) interactions, 2 for which the entropy is extensive even at zero temperature due to the exponential number of ground-state configurations.
Studying these systems at finite temperatures using Monte Carlo simulations typically poses no major challenges: An estimate of an observable (e.g., the magnetization) has to be measured and averaged over Monte Carlo time. This means that an average over different degenerate states is automatically taken into account. However, when studying the ground-state properties of such systems, it is imperative to ensure that an average over either all ground-state configurations or at least an unbiased subset is taken. While some algorithms are known to do a relatively fair sampling of the ground-state configurations (e.g., simulated annealing), 3,4 others, such as quantum annealing, [5][6][7][8] have been shown to favor certain configurations and exponentially suppress others.
Although finding ground-state energies is not typically harder with discrete energy distributions than in models with continuous energies, determining thermodynamic quantities requires thorough checks of the algorithms used. Numerically exact results are achievable only in limited scenarios, [9][10][11] though comparisons with theoretical predictions also may provide confidence that the technique is behaving correctly. 12,13 Here, we investigate the effectiveness of a heuristic optimization algorithm previously developed 14 for finding all ground states of a system, or, if that number is too large, finding a representative and unbiased sample, applied to a p-spin generalization of the d-dimensional Edwards-Anderson Ising spin-glass model. 15,16 The general applicability of this algorithm implies that it is useful for studying a wide range of hard problems where it is desirable to average over many degenerate opti-mal cases. We illustrate the usefulness of the method for studying spin-glass models by computing the zerotemperature ferromagnet-to-spin-glass phase transition as the disorder strength is varied in the disordered two-dimensional three-spin Baxter-Wu model and in the three-dimensional two-spin Edwards-Anderson model with bimodal interactions.
Using a moderate numerical effort, we show for the case of the three-dimensional two-spin Edwards-Anderson Ising spin-glass model that the critical threshold q c between a ferromagnetic and a spin-glass phase is q c = 0.2253 (7), in agreement with previous studies, 17 albeit with considerably smaller error bars.
In Sec. II we present an overview of the algorithm, followed by results on the two-dimensional three-spin model in Sec. III and results on the three-dimensional Edwards-Anderson spin glass in Sec. IV.

II. NUMERICAL METHOD
The Edwards Anderson model in three space dimensions has been shown to be NP-hard, 18 as has the spin glass with three-body interactions in two dimensions. 14 It is therefore expected that there are no exact algorithms capable of simulating large instances of either model at zero temperature. The genetic algorithm presented in Ref. 14 heuristically generates many optimal groundstate configurations of disordered p-spin models. The method starts from a population of N p random spin configurations and combines a highly-effective local search heuristic with triadic crossover genetic updates. 19 Although very large systems are unattainable with such an approach, it was shown to produce ground-state energies in the two-dimensional three-spin model with high confidence for N 42 2 = 1764 spins. 14 Genetic algorithm: The algorithm proceeds according to the following steps (a) Initialize the spins to a high-temperature configuration where all spins are assigned randomly.
(b) Perform a local search update to generate a low-energy configuration.
(c) If this configuration is already present in the population, go back to 1(a).
(d) Otherwise, accept this configuration and proceed to the next one.
2. Perform a triadic crossover update to generate two new configurations, C 1 and C 2 .
(a) Select three parents P 1 , P 2 and P 3 randomly.
(c) Select two configurations X and Y from the population at random.
3. Repeat step 2 N steps times or until all states of the system have the same energy.
Unlike in many genetic algorithms, mutations are not necessary to achieve good results with this method. Local search component: The local search algorithm used in step 1(b) above (described in detail in Ref. 14) generates clusters of nearby spins for which the energy is lowered if all spins in the cluster are flipped. This is achieved by performing a depth-first search starting from each spin in the system. The search algorithm takes two parameters, E max ≥ 0 and d max . At each step of the search, the spins adjacent to the current site that are not already in the cluster are considered to be added. As each spin is added, the energy of flipping the current spin and all its ancestors in the search tree (all current members of the cluster) is maintained. If the energy is negative, the algorithm has found an energy-lowering move, so all spins in the current cluster are flipped and the search terminates successfully. If the energy exceeds the threshold energy E max or the search depth becomes more than d max , this search branch is discarded. For E max and d max both large, this search algorithm gives strong results, but runs slowly. For use as part of the genetic algorithm above, speed is more important than finding the lowest energy, therefore we search for small clusters with few barriers, setting E max = 4 (the smallest possible change in energy in our model) and d max = L, where L is the linear size of the system. With N p and N steps sufficiently large, this genetic algorithm produces ground states with high confidence. 14 In general, to study a ferromagnet-paramagnet transition, it is useful to study the finite-size scaling of dimensionless quantities based on the magnetization of the system. By computing quantities such as the Binder ratio, 20 one can straightforwardly measure the location of the critical point, as well as the critical exponent ν. When the disorder distribution is continuous, each disorder sample has a unique ground state (up to global degeneracies), so finding this ground-state configuration immediately yields both the energy and magnetization; computing its energy and magnetization are therefore typically of comparable difficulty. But in the discrete-disorder case, the number of ground-state configurations N g is typically very large. Computing the magnetization requires either finding all ground states, or generating a uniform sampling to find a subset of these. This heuristic algorithm does not stop after it finds one ground-state configuration; it continues finding others until either all ground states are found or the algorithm can find no more ground states due to a built-in constraint.
For a genetic algorithm, the quality of the results depends both on the run time of the algorithm and the number of configurations being concurrently updated, the population size N p . For sufficient run time, our numerical results show that this algorithm finds all the ground-state configurations of a degenerate system if the number of ground states N g < N p . When N g > N p , which will inevitably be the case for some large systems, the probabilities for finding each ground state vary somewhat, although our tests suggest that the probability of finding a particular ground state is within some factor times the expected probability with a uniform distribution, with that factor typically being of order unity. For one instance of the 3-spin model studied in Sec. III with 600 ground states, we reran our algorithm 10 5 times to find the relative probability for the genetic algorithm to find each ground state. Note that disorder instances with this number of ground states are rare for L = 6; the median number of ground states we find over all disorder instances is 4. For this sample, when N p = 2N g all ground states were found in every run, and for N p = N g almost all ground states were found in nearly every runs. When N p = 0.48N g , the ground state with the largest basin of attraction is roughly eight times more likely to be seen than the ground state with the smallest. To probe the system size dependence, we also looked at a larger system of size L = 18, where we chose a sample that also has 600 ground states. In this case, the algorithm selects among those ground states much more uniformly, with only a factor of 2 between most and least likely ground states. For L = 18, the median number of ground states is 768; with N p > 5000, a large majority of the disorder instances considered have N g < N p and are expected to be treated fairly by the algorithm. These results are shown in Fig. 1.
Even when N g > N p , this performance strongly outperforms quantum annealing [5][6][7]21 where in the case of the fully-frustrated Villain model 22 certain ground-state configurations are exponentially suppressed. 8 We have : (Color online) Normalized relative probability p of each ground state occurring for varying population sizes. A single instance of disorder with N = 6 × 6 spins is considered, in addition to one instance of disorder with N = 18 × 18. The disorder cases were selected, for comparison, to have the same number of ground states, with Ng = 600. The ground states (indexed by i) are sorted according to their probability of occurring, giving a monotonically increasing plot. Even with Np < Ng, the probability for finding each ground state is a factor of order unity times the uniform probability, implying that one can do a reasonable approximation of the T = 0 thermodynamic state even when Np < Ng. With other factors being equal, the larger system size shows a more uniform sampling.
not found any trends in the magnetizations produced by the algorithm; using only this magnetization data, the ground states appear to be chosen in random order, even when N g > N p , as is shown in Sec. III, and Fig. 2. An average over the configurations that have been found therefore appears to be representative of the entire ensemble for this algorithm.

A. Model
The three-spin Ising model with random plaquette interactions is given by the Hamiltonian where the sum is over all triangular plaquettes in a twodimensional triangular lattice and each term involves the product of the Ising spins S i △ ∈ {±1} on the vertices of each plaquette. Without any loss of generality we set J ≡ 1. The plaquette couplings η △ are chosen independently and randomly according to a bimodal distribution, that is, with q = 0 being the pure ferromagnetic case. As q increases, the model goes through a phase transition from ferromagnet to paramagnet. A previous study at finite temperature using Monte Carlo methods finds this transition to occur at q c = 0.109 (2). 23 This model is of interest for a number of reasons. pspin models have been related to structural glasses. [24][25][26][27] Furthermore, when computing the error threshold of topological color codes 23,28 to bit-flip errors the problem maps 29 onto the 3-spin model in two space dimensions. Although the error threshold is given by the crossing of the ferromagnetic-paramagnetic phase boundary line with the Nishimori line, 30 the zero-temperature phase boundary delivers a lower bound for the threshold that is often easier to compute than to equilibrate finitetemperature Monte Carlo simulations.
The phase diagram for the Ising model on a square lattice shares much in common with the three-spin Ising model on a triangular lattice. The three-spin Ising model with uniform ferromagnetic interactions is exactly solvable and known as the Baxter-Wu model, 31 and the transition temperature on a triangular lattice is identical to that of the two-dimensional Ising model on a square lattice. Although the models are in different universality classes, 23,32 the two phase diagrams are quite similar even in the presence of disorder. Numerically, it appears that the multicritical points for these two models are very close to one another or identical. 23 The results presented here suggest that this similarity does not extend to the low-temperature regime below the Nishimori line, 30 with possible reentrant behavior being much weaker, if it exists at all, in this model.
In contrast to the Ising model with p = 2, which has a two-fold (global spin flip) degeneracy, this model is fourfold degenerate. The magnetization in the ordered state is not symmetric about zero, although the average over the four degeneracy sectors related by symmetry operations is zero: in the pure ferromagnet, one of the four degenerate ground states has M = 1, while the other three states each have M = −1/3. In the paramagnetic state there is no state which is much stronger than the others, and the disorder-averaged magnetization in each of these four sectors is zero. The zero-temperature thermodynamic state consists of all ground states with equal probability. To compute properties related to the magnetization, it is therefore not sufficient to find an arbitrary specific ground state of the model unless it can be demonstrated that it is chosen in an unbiased way.

B. Sampling results
When N p > N g , it is seen that the number of ground states found for the three-spin model is always a multiple of 4, and the average magnetization is 0. It is necessary that these both be true due to the fourfold degeneracy of the system. This result suggests that the algorithm is finding all ground states of the system. It is possible that there are other states which are not found, but we emphasize that the algorithm is given no information about the symmetries of the problem, so states which are related to one another by global spin-flip symmetries are not necessarily easy for the algorithm to link together. Thus is appears that we are seeing all ground states when N p > N g .
When the number of ground states is large, so that N g > N p , it also appears that the subset of all ground states that the algorithm finds is representative of the entire ensemble. To test this, we study the magnetization of the ground states. The results for two typical cases are shown in Fig. 2, where the magnetization is shown for the first 5000 ground states found in the order i that they are found by the algorithm. There are no apparent temporal correlations in the output and the running average quickly approaches zero, although there are some statistical fluctuations about this value because not all ground states have been found. This, in turn, justifies the use of this algorithm for sampling.
The two samples shown in Fig. 2 are for two different instances of disorder with the same disorder strength q = 0.105 and the same system size (N = 18 2 ). This disorder strength is near the ferromagnet to spin-glass phase transition but slightly below it, so that the system is still ferromagnetic. In the sample shown in the top panel, the states all resemble the pure three-spin spin-glass problem, with 1/4 of the configurations having M ≈ 1, and the other 3/4 having M ≈ −1/3. Unlike in the pure case, there are many distinct ground states, but each one follows the expectations for a ferromagnetic system. The sample depicted in the bottom panel shows the onset of disorder: while the majority of the states still have M ≈ 1 or M ≈ −1/3, a sizable fraction of the states here have M slightly below 0, and still others are between 0 and 1/2. This implies that there are large collective excitations with zero-energy, which is a signature of glassy/frustrated systems. Thus, even for the same disorder strength, the system will behave ferromagnetically for some instances of the disorder and appear to be more glassy for others.
It is instructive to compare against the twodimensional Ising spin-glass model with p = 2. Even though efficient algorithms have long been known for computing ground states in this case, 18 these techniques have not been directly useful in finding all of the exponentially-many ground states, or even in finding unbiased randomly chosen ground states. Heuristic measures have therefore been employed to find the groundstate magnetization, 33 although an efficient algorithm for fair sampling of the ground states for this model has been developed recently. 34 This model has been shown to exhibit reentrance from both zero-and finite-temperature simulations, 33 on the three-spin model below the Nishimori line for comparison.

C. Ferromagnet-to-spin-glass transition
Although this algorithm does not permit the computation of ground states with high confidence in very large systems, the intermediate system sizes accessible are adequate to precisely determine the location of the T = 0 ferromagnet-to-spin-glass phase transition as the disorder strength q varies. To do this, we compute the Binder ratio 20 for a restricted set of magnetizations. Due to the details of the four-fold degeneracy in this model including the lack of a global spin-flip symmetry, the ferromagnetic state consists of three states with weak negative magnetization for every highly magnetized state with positive magnetization. The paramagnetic state has no sector with stronger magnetization than the others. The fluctuations are strongest in the ferromagnetic sector, because there is less change from the behaviors in the other sectors between the two phases, so it is natural to focus on only the 1/4 of states that have a large magnetization. In Fig. 2, these are the states with M near 1. Let and for numbered colors k = 1, 2, 3 of the tripartite lattice, with r ik = 1 if site i has color k and −1 otherwise. We focus on the state of the system with the largest magnetization, calling this the restricted magnetization For any given spin configuration, it is straightforward to compute all four of these magnetizations. This is analogous to taking the absolute value for the p = 2 Ising model, such that only configurations in the same state are compared against one another. We then compute the restricted Binder ratio g R defined by where [· · · ] av is an average over samples and · · · represents an average over ground-state magnetizations for the sample. This is the same as the standard definition of the Binder ratio except that the restricted magnetization M R is used in place of M . For the p = 2 Ising model, taking the absolute value would not make any difference, but here the symmetric states have different magnetizations. The restricted Binder ratio is dimensionless and is 1 in the ferromagnetic phase and 0 in the paramagnetic phase (here, the T = 0 spin glass). In a finite system of linear scale L with disorder strength in the vicinity of q c , g R is expected to scale as Equation (8) allows for the extraction of the critical disorder strength q c as well as the correlation-length exponent ν. To estimate the optimal values of the critical parameters we perform a finite-size scaling analysis where we tune q c and ν until the chi-squared of a fit to a thirdorder polynomial in the vicinity of L 1/ν (q − q c ) 1 is minimized. Error bars are determined by a bootstrap analysis. 37,38 Our best estimates of the critical parameters are q c = 0.1072(1) ν = 1.5(1) .
The scaling collapse using these critical parameters is shown in Fig. 3. To verify our results, we have also computed the standard Binder ratio on the unrestricted magnetizations from the same data, which produced consistent results for both q c and ν. The only noticeable difference in the results is the form of the scaling function, that is, g R (q = q c ) ≈ 0.955, while g(q = q c ) ≈ 0.23. A summary of the simulation parameters is shown in Table I. Finally, to treat the possibility of finite-size corrections to scaling, we separate the components of the fit in Fig. 3 into L-2L pairs. We perform the curve-fitting procedure to extract q c using each pair of system sizes. The triangular lattice is tripartite only for L x mod 3 = 0 and L y mod 2 = 0, so these comparisons are only possible for and N = L×L triangular lattice if L is a multiple of 6. To increase the number of cases we can handle, we consider L x to be any multiple of 3, and if L x is not divisible by 2, we simulate for both L x × (L x − 1) and L x × (L x + 1) systems. This gives two results which are presumably biased in opposite directions, which can be seen for the odd system sizes in Fig. 3. However, our statistical errors are as large as the bias introduced, so we simply average the two results. In comparing sequences of L-2L pairs, we check if the result depends on L to extrapolate to the thermodynamic limit. We find that the extracted values of q c (L, 2L) show no visible dependence on system size, as shown in Fig. 4. Thus the above estimates of q c and ν appear to be indicative of the values in the thermodynamic limit. In contrast to the two-dimensional Ising model, the three-spin model appears to have much weaker reentrance, if it is reentrant at all. The phase diagrams for the two models therefore differ considerably below the Nishimori line.

IV. THREE-DIMENSIONAL ISING SPIN GLASS
As another application of this algorithm, we investigate the ferromagnet-to-spin-glass transition in the three-dimensional Edwards Anderson Ising spin glass with bimodal disorder. This model is important in the study of disordered materials. 2 It is also notoriously difficult to study numerically. It is NP-hard, 18 so exact methods are quite limited in scope, 39 and even highlysophisticated heuristic methods can only produce reliable results up to L = 14. 17

A. Model
The two-spin Edwards-Anderson Ising spin glass is given by the Hamiltonian where the spins S i sit at the sites of a cubic lattice and pairwise interactions (p = 2) are over nearest neighbors. Here we consider bimodal disorder, where the bond strengths are given by the same probability distribution as the plaquettes in Eq. (2), The zero-temperature ferromagnet-to-spin-glass transition has been investigated for this problem using a sophisticated heuristic algorithm, 17 where the critical disorder strength was found to be q c = 0.222 (5). This work did not take into account possible biases in the sampling of ground states, which could lead to incorrect results, as pointed out in Ref. 40. The finite-temperature multicritical point has been evaluated for this model in a number of studies, with the results summarized in Table II and shown graphically in Fig. 5. Note that, where necessary, we have converted estimates so they are all in the same terms. In particular, Ref. 41 quotes only the temperature T * c of the multicritical point, which is easily converted to q * c , and the value given for ν in Ref. 42 is for thermal perturbations, with ν = 1/y 2 , while the value we quote (1/y 1 ) is more relevant to low-temperature studies because it describes the disorder perturbations that are used at the zero-temperature fixed point. The temperature of the multicritical point (MCP) can be determined directly from the Nishimori condition 1 − 2q = tanh(1/T ), (12) with values in the literature in the range T * c ≈ 1.67-1.69. It is believed that the data below the multicritical point are universal and in a different universality class than at the multicritical point. 43 This view is consistent with the known data. Furthermore, the q c values show strong evidence for reentrance in the phase diagram.

B. Results
We use the same code for simulating the three-spin Ising spin glass to also find ground states of the p = 2 Ising spin glass in three dimensions. We have computed the Binder ratio in the vicinity of the critical disorder strength for systems of sizes 4 3 -10 3 . In the Edwards-Anderson model, the restricted Binder ratio is identical   Table II (tagged via roman numerals). The center of an ellipse corresponds to the optimal estimate of qc and ν, while the size of the ellipse corresponds to the error bar attached to it. The red ellipses [I -III] are at the multicritical point, which is believed be in a different universality class than the low-and zero-temperature case. The blue ellipses [VI, VII] are zero-temperature estimates (the filled ellipse represents the results from this study) and the purple ellipses [IV, V] are for finite temperatures below the multicritical point. The fact that the finite-temperature (red) and zero-temperature (blue) estimates do not overlap in the horizontal direction is indicative of a reentrant behavior in the phase diagram of the model.
to the standard definition of the Binder ratio because only even moments of the magnetization are used; here g = g R , so the scaling relation in Eq. (8) is again expected to hold. Figure 6 shows the scaling collapse for this model. The optimal values of the critical parameters  extracted from this best fit are q c = 0.2253(7) ν = 1.07 (7), (13) in agreement with the results of Ref. 17. The system sizes available do not permit the more thorough finite-size scaling analysis shown in Fig. 4, although the corrections to scaling appear to be quite weak, suggesting that this result is robust. A summary of the simulation parameters is shown in Table III.

V. CONCLUSIONS
We have demonstrated the applicability of the genetic ground-state algorithm presented in Ref. 14 to sample among the many ground states in highly-degenerate NPhard models. Using this technique, spin configurations and therefore magnetizations of typical ground-state configurations may be accessed, allowing the construction of parameters such as the Binder ratio without bias, that are effective for determining the location of phase transitions.
Above the Nishimori line, the phase diagram of the three-spin Ising model on a triangular lattice closely resembles that of the square-lattice Ising model with twospin interactions. Although the models are in different universality classes, the measured critical values for the pure case and on the Nishimori line are the same. The results here show that this similarity does not extend below the Nishimori line, that is, the models have quite different transition disorder strengths at T = 0 (see Fig. 4). We show that q c (T = 0) = 0.1072(1).
Furthermore, using a modest numerical effort we determine the ferromagnet-to-spin-glass transition for the three-dimensional Edwards-Anderson Ising spin glass, that is, q c = 0.2253 (7), in agreement with previous measurements. 17 Comparing against estimates of the multicritical point, this three-dimensional model pos-sesses reentrance in its phase diagram by an amount which is similar to the two-dimensional model (similar results were obtained recently in Ref. 43).