A dedicated algorithm for calculating ground states for the triangular random bond Ising model
aa r X i v : . [ c ond - m a t . d i s - nn ] J u l A dedicated algorithm for calculating ground states for the triangular random bond Isingmodel
O. Melchert, A. K. Hartmann
Institut f¨ur Physik, Universit¨at Oldenburg, Carl-von-Ossietzky Straße 9–11, 26111 Oldenburg, Germany
Abstract
In the presented article we present an algorithm for the computation of ground state spin configurations for the 2 d random bondIsing model on planar triangular lattice graphs. Therefore, it is explained how the respective ground state problem can be mappedto an auxiliary minimum-weight perfect matching problem, solvable in polynomial time. Consequently, the ground state propertiesas well as minimum-energy domain wall (MEDW) excitations for very large 2 d systems, e.g. lattice graphs with up to N = × T = Keywords:
Random bond Ising model, Negative-weight percolation, Groundstate phase transitions
1. Introduction
Triggered by the exchange of ideas between computer sci-ence and theoretical physics in the past decades, it was realizedthat several basic problems in the context of disordered systemsrelate to “easy” optimization problems. These are problemswhere the solution time is polynomial in the size of the problemdescription. As a result, many disordered systems can now beanalyzed numerically exact through computer simulations byusing fast combinatorial optimization algorithms [1, 2, 3]. E.g.,ground state (GS) spin configurations for the random-field Isingmagnet (in any dimension d ) can be obtained by computing themaximum flow for an auxiliary network problem [2]. Anotherexample is the 2 d Ising spin glass (ISG), where the lattice canbe embedded in a plane. For this model, the problem of find-ing a GS spin configuration for a given realization of the nearestneighbor couplings can be mapped to an appropriate minimum-weight perfect-matching (MWPM) problem [2, 4]. Finally, theMWPM problem can be solved in polynomial time by meansof exact combinatorial optimization algorithms [5]. Thus, theplanar 2 d ISG can be studied directly at zero temperature with-out equilibration problems and within polynomial time. Hence,very large systems can be considered, giving very precise andreliable estimates for the observables. Actually, there are dif-ferent approaches that allow for an exact computation of GSsfor the planar 2 d ISG [6, 7, 8]. Albeit all of these approachesrely on the computation of MWPMs on an auxiliary graph, theydi ff er regarding the subtleties of the mapping to the respectiveauxiliary problem. The most e ffi cient of these approaches (see Email addresses: [email protected] (O.Melchert), [email protected] (A. K. Hartmann)
Ref. [8]) is based on the Kasteleyn treatment of the Ising model[9], which previously was also used to obtain extended groundstates for the 2 d ISG with fully periodic boundary conditions[10].Here, we introduce a dedicated algorithm that yields exactGS spin configurations for the 2 d random-bond Ising model(RBIM) on planar triangular lattice graphs. As the previousapproaches, the algorithm presented here requires to solve anassociated MWPM problem. The corresponding mapping usesa relation between perfect matchings and paths on a graph[11, 12]. In e ff ect, these paths can be used to partition thegraph into domains of up and down spins that comprise a GSspin configuration, see Fig. 1. Consequently, the GS propertiesas well as minimum-energy domain wall (MEDW) excitations,see Fig. 1, can be analyzed very fast. The presented algorithmenables us to study large systems, while allowing for an appro-priate disorder average within a reasonable amount of comput-ing time. In this regard, it requires to compute a MWPM for anauxiliary graph with O ( N ) edges only (wherein N is the num-ber of spins on the lattice). However, note that the algorithmpresented here is asymptotically not faster than the algorithmpresented in Refs. [8, 10], but it highlights the algorithmic rela-tion between the GS problem for spin glasses and the recentlyproposed negative-weight percolation (NWP) problem [13].In the presented article, we investigate the critical behaviorof the T = d RBIM, signaled by a breakdown of the magnetization,using finite-size scaling (FSS) analyses of the MEDW excita-tion energy. In this regard, we obtain a highly precise estimateof the critical point for the triangular lattice geometry and weverify the critical exponents obtained earlier for the RBIM onthe planar square lattice [14, 15]. Finally, we contrast our nu-
Preprint submitted to Elsevier June 7, 2018 b)(a)(c)
Figure 1: Samples of a ± J random bond Ising spin system on a triangularlattice of side length L =
32. The samples are taken at three di ff erent values ofthe disorder parameter p . (a) p = .
15 characterized by a ferromagnetic GS, (b) p = . p = .
5, i.e. the canonical ± J ISG. In the figure, periodic BCs are indicated by the dashed vertical lines.From left to right: (Left) Transition graph that describes the di ff erence betweena ferromagnetic reference spin configuration and the GS spin configuration,(center) corresponding GS, (right) MEDW excitation relative to the GS. merical results with previous simulations and presumably exactresults [16].The remainder of the presented article is organized as fol-lows. In section 2, we introduce the model in more detail andwe outline the algorithm used to compute the GS spin config-urations. In section 3, we present the results of our numericalsimulations and in section 4 we conclude with a summary.
2. Model and Algorithm
In the presented article we perform GS calculations for the 2 d RBIM, where the respective model consists of N = L × L Isingspins σ = ( σ , . . . , σ N ), where σ i = ±
1, located on the sitesof a planar triangular lattice graph. Therein, the energy of agiven spin configuration is measured by the Edwards-AndersonHamiltonian H ( σ ) = − X h i , j i J i j σ i σ j , (1)where the sum runs over all pairs of nearest-neighbor spins (onthe triangular lattice) with periodic boundary conditions (BCs)in the x -direction and free BCs in the y -direction. In the aboveenergy function, the bonds J i j are quenched random variablesdrawn from the disorder distribution P ( J ) = p δ ( J + + (1 − p ) δ ( J − . (2)Therein, one realization of the disorder consists of a randomfraction p of antiferromagnetic bonds ( J = −
1) that prefer anantiparallel alignment of the coupled spins, and a fraction (1 − p ) of ferromagnetic bonds ( J =
1) in favor of parallel alignedspins. In general, the competitive nature of these interactionsgives rise to frustration. A plaquette, i.e. an elementary trian-gle on the lattice, is said to be frustrated if it is bordered by anodd number of antiferromagnetic bonds. In e ff ect, frustrationrules out a GS (i.e. a minimizer σ GS of Eq. 1) in which all thebonds are satisfied. As limiting cases one can identify the Isingferromagnet at p = ± J ISG at p = / p we expect tofind a ferromagnetic phase (spin-glass phase) for p < p c ( p > p c ),wherein p c denotes the critical point at which the T = N [17, 18]. Apart fromthe GSs, we here also aim to characterize the energetic proper-ties of MEDW excitations. A domain wall is an interface thatspans the system in the direction with thee free BCs. Now, theMEDW is such an interface with an excitation energy δ E thatis minimal among all possible domain walls. Due to the exten-sive degeneracy of the GSs, the “lattice-path” associated with aMEDW is not unique. I.e., there are many DWs with minimalexcitation energy. Albeit the geometric properties of a MEDWare not unique, its excitation energy is unique. MEDWs forthree di ff erent values of the disorder parameter p are illustratedin Fig. 1.We now give a brief description of the algorithm that we useto compute the GSs. Therefore, we first set a reference spinconfiguration σ R . The most convenient choice is a maximallypolarized, i.e. ferromagnetic, configuration σ R = ( + , . . . , + L extra nodes ontop and at the bottom of the dual in order to account for thefree BCs along that direction and to maintain the honeycombstructure of the respective graph. This means, a triangular spinlattice of size L × L is transformed to a honeycomb lattice withan over all number of (2 L ) × ( L +
1) nodes. Hence, the topolog-ical dual graph associated to the triangular spin lattice is modi-fied to some extend. Further note that adjacent extra nodes areconnected by edges e that carry a weight ω ( e ) =
0. All otheredges e on the dual graph get an edge-weight ω ( e ) ≡ J i j σ R , i σ R , j .Therein, e is assumed to cross a bond J i j on the spin lattice,where J i j couples the two spins σ R , i / j , see Fig. 2(a). Conse-quently, the edge weight on the weighted dual is positive (neg-ative), if the corresponding bond on the spin lattice is satisfied(broken) with respect to σ R . A pivotal observation is, that thereexists an equivalence between clusters of adjacent spins on thespin lattice that might be flipped in order to decrease the config-urational energy of σ R and negative-weighted loops (i.e. closedpaths) on the weighted dual graph. In this regard, if a loopwith negative weight on the dual is found, the cluster of spinssurrounded by this loop can be flipped so as to decrease theconfigurational energy of σ R . Finally, to obtain a GS spin con-figuration one needs to find a minimum-weight set of negative-weighted loops on the dual graph (see discussion below). Thisset of loops comprises the transition graph for the given real-2 a) (b) (c) w(e)=1w(e)=0J =1 ij Figure 2: Illustration of the computation of the transition graph that allows to determine a GS spin configuration on a planar triangular RBIM. The figure illustratesa sample system of side length L = G (grey edges, triangular geometry), wherea solid (dashed) line indicates a ferromagnetic (antiferromagnetic) bond, to the weighted dual graph G D (black edges, honeycomb geometry). Note that additionaledges (dotted lines) where introduced to account for the open boundary conditions in the respective direction. These edges carry zero weight. (b) Minimum weightset of loops (bold black edges) on G D as obtained from a mapping to the NWP problem (not shown here, see Ref. [13]). (c) Loop on the dual surrounding a cluster ofspins on the spin lattice. If the orientation of the spins is chosen as explained in the text, this procedure yields a GS spin configuration. In the figure, spin orientationsare distinguished by gray filled and non-filled circles. ization of the disorder, as illustrated in Fig. 2(b). Since initiallya ferromagnetic reference configuration was chosen, the GS isobtained if the orientation of the spins on the spin lattice is cho-sen such, that (i) spins within a cluster are aligned in the samedirection and (ii) spins in adjacent clusters are aligned in oppo-site directions. The resulting GS is indicated in Fig. 2(c), seealso Fig. 1.Here, for the 2 d RBIM on a planar triangular lattice, wherethe dual has a honeycomb geometry, the minimum weight setof loops on the weighted dual can be obtained by means of amapping to the NWP problem, as explained in [13]. In brief,the NWP statement consists in the task to find a minimum-weight set of nonintersecting negative-weighted loops for agiven weighted graph. Therefore, it considers a minimum-weight perfect matching problem on an associated auxiliarygraph with O ( N ) edges (provided that the input graph has O ( N )edges, as it is the case here), from which the set of loops canbe deduced. Note that the mapping to the NWP problem yieldsthe correct transition graph only for this particular lattice setup,since any two-coloring of the spin lattice (i.e. assignment ofup / down spin orientations) can be composed by loops on thedual that do not intersect. That means, each site on the dualis an end-node of either 0 or 2 loop segments, as e.g. in Fig.2(b). In contrast, two-colorings of the spins on a square lat-tice might involve loops on the dual that involve figure-8 twists.That means, each site on the dual is end-node of either 0, 2 or 4loop segments. For the latter problem, a di ff erent mapping [8]was used recently to obtain exact GSs for 2 d ISGs on a squarelattice with free BCs in at least one direction within polynomialtime. This mapping was further used to compute “extended”GSs for the 2 d ISG with fully periodic BCs [10].Now, the interpretation of the T = p , there are only few bonds on the spin lattice that are notsatisfied by the reference spin configuration. Accordingly, thereare only few small loops that comprise the transition graph. Forall nonzero values of p , a su ffi ciently large lattice will featureat least some small loops that surround an elementary plaquette on the dual. These small loops correspond to local “manipu-lations” of the order parameter (i.e. the magnetization), only.Hence, in the thermodynamic limit, the GS has still ferromag-netic order (see Fig. 1(a)). However, if the value of p increasesand exceeds a critical value p c , large loops appear that have alinear extension of the order of the system size and eventuallyspan the system along the direction with the periodic boundaryconditions. These loops represent global manipulations of theorder parameter, that, in the thermodynamic limit, destroy theferromagnetic order of the GS (see Figs. 1(b),(c)).Once we obtained a GS spin configuration in this manner, wecompute a MEDW by means of a similar mapping, thoroughlyexplained in Ref. [19]. In the following we will use the proce-dure outlined above to obtain GSs and to investigate MEDWsfor the RBIM introduced above.
3. Results
As pointed out above, at small values of p there exists anordered ferromagnetic phase, while for large values of p a spin-glass ordered phase appears. A proper order parameter to char-acterize the respective FM-SG transition is the magnetizationper spin m L = | P i σ i | / L for a system of linear extend L . Be-low, we perform a finite-size scaling analysis (FSS) in order tolocate the critical point p c and to estimate the critical exponentsthat describe the scaling behavior of the magnetization in the Table 1: Critical exponents for the 2 d RBIM. From left to right: Problem setup(SQ = square lattice, TR = triangular lattice), critical exponent of the correlationlength ν , order parameter exponent β , and exponents φ and ψ that characterizethe scaling of the MEDW excitation energy. The figures for SQ-a are taken fromRef. [15]. The figures for SQ-b are taken from Ref. [14]. Setup ν β φ ψ SQ-a 1 . . . . . . . . . . . . b L ( p ) p ISG-tr- ± JL=2432486496128 0.4 0.6 1-0.4 -0.2 0 b L ( p ) (p-p c )L ν p c =0.1584(3) ν =1.47(6) Figure 3: Results of the finite-size scaling analysis for the binder parameter b L ( p ), considering di ff erent system sizes L . The main plot shows the unscaleddata close to the critical point and the insets illustrate the data collapse obtainedafter rescaling the raw data using the scaling assumption discussed in the textand scaling parameters as listed in Tab. 1. vicinity of the critical point. Therefore, we first consider theBinder parameter [20] b L = (cid:18) − h m L ih m L i (cid:19) (3)associated with the magnetization. It is expected to scaleas b L ( p ) ∼ f [( p − p c ) L /ν ], wherein f [ · ] signifies a size-independent scaling function and ν denotes the critical expo-nent that describes the divergence of the correlation length asthe critical point is approached. Here, we simulated triangularsystems of side length L =
24 through 128 at various valuesof the disorder parameter p . Observables are averaged over64 000 samples for the largest systems and we used the datacollapse generated by the scaling assumption above to obtain p c = . ν = . S = .
94 of thedata collapse [21, 22], see Fig. 3. In general, the above scalingrelation holds best near the critical point and one can expect thatthere are corrections to scaling o ff criticality. As a remedy, werestricted the latter scaling analysis to the interval [ − . , + . h m L ( p ) i ∼ L − β/ν f [( p − p c ) L /ν ], where f [ · ] denotes a size-independent function, andwhere the order parameter exponent β can be obtained after fix-ing ν and p c to the values obtained from the analysis of theBinder parameter. The best data collapse ( S = .
01) was ob-tained for the choice β = . h δ E i according to the scaling assumption h δ E i ∼ L ψ f [( p − p c ) L φ ], see Ref. [23], yields the critical point p c = . -1 -0.4 -0.2 0 0.2 L − ψ 〈 δ E 〉 (p-p c )L φ ISG-tr- ± JL=24324864961280.1110 0.1 0.15 0.2 〈 δ E 〉 p Figure 4: Results of the finite-size scaling analysis for the average MEDWenergy h δ E i . The main plot shows the data collapse after rescaling the raw datausing the scaling parameters listed in Tab. 1 and the inset illustrates the unscaleddata close to the critical point. parameter. The critical exponents ψ and φ are listed in Tab.1. Therein, we restricted the scaling analysis to the interval[ − . , + .
1] and obtained a best data collapse with S = . p c we performed additional simu-lations for spin lattices of up to 384 ×
384 spins (and 3 . × samples), i.e. weighted dual graphs of up to 768 ×
385 nodes.Upon analysis of the data we obtain the estimate β = . h m i ∼ ( L + ∆ L ) − β/ν ,wherein ∆ L = O (1). Considering the scaling of the averageMEDW excitation energy h δ E i and using a similar scaling as-sumption as above, we found ψ = . ff er regarding theirgeometric properties. However, here we also analyze the av-erage length h ℓ i of the particular MEDWs obtained within thesimulations, see Fig. 5(b). Therefore, we considered a scal-ing according to the form h ℓ i ∼ ( L + ∆ L ) d f , wherein d f signi-fies the fractal dimension of the MEDWs at p c . We obtained d f = . ∆ L = O (1)), which is in agreement withthe value d f = . T = d square lattice, see Ref. [15].
4. Conclusions
In the presented article we have illustrated how GSs for the2 d RBIM on planar triangular lattice graphs can be computedby a mapping to the NWP problem. I.e., the problem of find-ing a GS spin configuration for a planar 2 d triangular RBIM isequivalent to the NWP problem on a properly weighted corre-4 L 〈 m 〉 (a)
24 96 384 L 〈 ℓ 〉〈δ E 〉 (b) Figure 5: Results of the finite-size scaling analysis at the critical point, wherethe T = h m i , and, (b) scaling of the average MEDW length h ℓ i and the MEDW excita-tion energy h δ E i . sponding dual graph that exhibits a honeycomb structure. Us-ing this approach, we have investigated GSs and MEDW exci-tations for the respective lattice structure. Therein, a disorderparameter could be used to distinguish a ferromagnetic and aspin-glass ordered phase. We characterized the corresponding T = d RBIM on a planar square lattice by consider-ing a Gaussian bond distribution with ferromagnetic bias [15]or a bimodal bond distribution [14], as listed in Tab. 1.Hence, the results for the triangular lattice structure obtainedhere highlights the universality of the T = p c and ν found here agree well with the values p c = . ν = . d lattice graphs with a hon-eycomb geometry and fully periodic boundary conditions [13].Finally, the location of the critical point obtained here via FSSanalysis is close to the theoretical prediction p c , tr = .
15, thatwas obtained for systems with fully periodic boundary condi-tions using the adjoined problem approach [24, 16].
Acknowledgment
We are grateful to B. Ahrens for a critical reading of themanuscript. OM acknowledges financial support from the Volk-swagenStiftung (Germany) within the program “Nachwuchs-gruppen an Universit¨aten”. The simulations were performed atthe GOLEM I cluster for scientific computing at the Universityof Oldenburg (Germany).
References [1] S. Bastea, A. Burkov, C. Moukarzel, P. M. Duxbury, Combinatorial op-timization methods in disordered systems, Comp. Phys. Comm. 121- 122 (1999) 199. Proceedings of the Europhysics Conference on Compu-tational Physics CCP 1998.[2] A. K. Hartmann, H. Rieger, Optimization Algorithms in Physics, Wiley-VCH, Weinheim, 2001.[3] H. Rieger, Application of exact combinatorial optimization algorithms tothe physics of disordered systems, Comp. Phys. Comm. 147 (2002) 702.[4] A. K. Hartmann, Domain walls, droplets and barriers in two-dimensionalising spin glasses, in: J. W (Ed.), Rugged Free Energy Landscapes,Springer, Berlin, 2007, pp. 67 – 106.[5] W. Cook, A. Rohe, Computing minimum-weight perfect matchings, IN-FORMS J. Computing 11 (1999) 138–148.[6] I. Bieche, R. Maynard, R. Rammal, J. P. Uhry, On the ground states of thefrustration model of a spin glass by a matching method of graph theory,J. Phys. A: Math. Gen. 13 (1980) 2553–2576.[7] F. Barahona, On the computational complexity of Ising spin glass models,J. Phys. A: Math. Gen. 15 (1982) 3241–3253.[8] G. Pardella, F. Liers, Exact ground states of large two-dimensional planarIsing spin glasses, Phys. Rev. E 78 (2008) 056705.[9] P. W. Kasteleyn, Dimer Statistics and Phase Transitions, J. Math. Phys. 4(1963) 287–293.[10] C. K. Thomas, A. A. Middleton, Matching Kasteleyn cities for spin glassground states, Phys. Rev. B 76 (2007) 220406.[11] R. K. Ahuja, T. L. Magnanti, J. B. Orlin, Network Flows: Theory, Algo-rithms, and Applications, Prentice Hall, 1993.[12] O. Melchert, PhD thesis, not published, 2009.[13] O. Melchert, A. K. Hartmann, Negative-weight percolation, New. J. Phys.10 (2008) 043039.[14] C. Amoruso, A. K. Hartmann, Domain-wall energies and magnetizationof the two-dimensional random-bond Ising model, Phys. Rev. B 70 (2004)134425. Note that the value of β = . β = . T = ± J Isingmodels, Physica A 245 (1997) 560 – 574.[17] L. Saul, M. Kardar, The 2 d ± J Ising spin glass: exact partition functionsin polynomial time, Nucl. Phys. B 432 (1994) 641.[18] J. W. Landry, S. N. Coppersmith, Ground states of two-dimensional ± J Edwards-Anderson spin glasses, Phys. Rev. B 65 (2002) 134404.[19] O. Melchert, A. K. Hartmann, Fractal dimension of domain walls in two-dimensional Ising spin glasses, Phys. Rev. B 76 (2007) 174411.[20] K. Binder, Finite size scaling analysis of ising model block distributionfunctions, Z. Phys. B 43 (1981) 119.[21] J. Houdayer, A. K. Hartmann, Low-temperature behavior of two-dimensional Gaussian Ising spin glasses, Phys. Rev. B
70 (2004) 014418. S measures the mean square distance of the scaled data to the mastercurve in units of standart errors.[22] O. Melchert, autoScale.py - A program for automatic finite-size scalinganalyses: A user’s guide, Preprint: arXiv:0910.5403v1 (2009).[23] N. Kawashima, H. Rieger, Finite-size scaling analysis of exact groundstates for ± J spin glass models in two dimensions, Europhys. Lett. 39(1997) 85–90.[24] J. Bendisch, Groundstate threshold pc in ising frustration systems on 2dregular lattices, Physica A 202 (1994) 48 – 67.spin glass models in two dimensions, Europhys. Lett. 39(1997) 85–90.[24] J. Bendisch, Groundstate threshold pc in ising frustration systems on 2dregular lattices, Physica A 202 (1994) 48 – 67.