Negative-weight percolation

We describe a percolation problem on lattices (graphs, networks), with edge weights drawn from disorder distributions that allow for weights (or distances) of either sign, i.e. including negative weights. We are interested whether there are spanning paths or loops of total negative weight. This kind of percolation problem is fundamentally different from conventional percolation problems, e.g. it does not exhibit transitivity, hence no simple definition of clusters, and several spanning paths/loops might coexist in the percolation regime at the same time. Furthermore, to study this percolation problem numerically, one has to perform a non-trivial transformation of the original graph and apply sophisticated matching algorithms. Using this approach, we study the corresponding percolation transitions on large square, hexagonal and cubic lattices for two types of disorder distributions and determine the critical exponents. The results show that negative-weight percolation is in a different universality class compared to conventional bond/site percolation. On the other hand, negative-weight percolation seems to be related to the ferromagnet/spin-glass transition of random-bond Ising systems, at least in two dimensions.


Introduction
Percolation [13] is one of the most-fundamental problems in statistical mechanics. Many phase transitions in physical systems can be explained in terms of a percolation transition. The pivotal property of percolation is connectivity. One can describe this in terms of weighted graphs (often also called networks) with nonzero or zero weights, corresponding to occupied/connected or unoccupied/disconnected edges. Here, we extend the problem to negative weights, in particular we study bond percolation. As an example, one can imagine an agent traveling on a graph, who has, when traversing an edge, either to pay some resource (positive weight) or he is able (once) to harvest some resource (negative weight). Paths including negative edge weights also appear in the context of domain walls in random-bond Ising systems [10]. One percolation problem is whether there exists a path (or loop) spanning the full system with negative total weight, such that each edge is traversed at most once. Hence, the percolating objects are paths and it is sufficient to look for minimum-weighted (or "shortest") paths (MWPs). Percolation properties of string-like objects have been studied occasionally [7,12,5,14,4,11], but to our knowledge never allowing for negative weights. One realizes immediately that the negative weights lead to properties which are fundamentally different from conventional percolation. E.g., negative-weight percolation (NWP) lacks transitivity: If there is a valid (negativeweight) path S → A → B from S to B via A, it is possible that the path S → A is not valid, as for the paths 0−1−2−3/0−1 in figure 1. Hence, there is no simple definition of percolation clusters in this case. Note that, since strings are thus considered, several percolating strings or loops might coexist in the same sample. This we observe indeed, see below.
Here, we study NWP numerically. First we introduce the model. Then we outline the algorithm to obtain MWPs, since no standard shortest-path algorithm can be used. The results for different types of two-(2d) and three-dimensional (3d) systems, in particular the critical exponents describing the percolation transitions, show that NWP is fundamentally different compared to conventional percolation and from previous models with string-like percolation. Thus, the study of percolation models exhibiting negative weights might lead to many new insights concerning the behavior of disordered systems.

Model and Algorithm
We consider 2d and 3d lattice graphs G = (V, E) with side length L, where adjacent sites i, j ∈ V are joined by undirected edges e = (ij) ∈ E. Weights ω(e) are associated with the edges, representing quenched random variables. We consider either bimodal (±J) or "Gaussian-like" (Gl) distribution of the edge weights, where ρ denotes the fraction of negative or Gaussian-distributed edge weights, respectively, among edges with unit weight (fraction 1 − ρ). The bimodal weights are taken to be ±1 and the Gaussian weights have zero mean and unit width. These weight distributions explicitly allow for loops L with negative weight, given by ω L ≡ e∈L ω(e). For any non-zero value of ρ, a sufficiently large systems will exhibit at least small loops with negative total weight. We investigate (i) MWPs in the presence of negative weighted loops on square lattices with periodic boundary conditions (BCs) in one and free BCs in the other direction. The path ends are allowed to terminate on the free boundaries. Further, we study (ii) minimum-weighted configurations of negative-weighted loops on square, 2d hexagonal and cubic lattices with fully periodic BCs, see figure 2. In either case we find critical values of ρ above which (i) negative-weight paths appear that span the lattice across the direction with periodic BCs or (ii) percolating loops emerge, that span the lattice along some direction.
Since we are looking for MWPs on undirected graphs in the presence of negative weights, traditional "shortest path" algorithms cannot be applied. The reason is that for applying these traditional algorithms a special condition must hold: For the distance d(i) of any shortest path from a source 0 to i = 0, d(i) = min j∈N (i) (d(j) + ω(j, i)) holds, where N (i) denotes the set of neighbors of i. This equation is not fulfilled in our case, as can be seen from figure 1 for node 4. It has minimum distance −1 to node 0, but it is connected to node 5 via an edge of weight ω(5, 4) = −1 and d(5) = −2. Hence, a different approach has to be applied. We determine MWPs and loop configurations through an appropriate transformation of the original graph, detailed in [2], and obtaining a minimum-weighted perfect matching (MWPM) [6,8], by using exact combinatorial optimization algorithms.
Here we give a concise description of the mapping, pictured as 3-step procedure illustrated in figure 3: (1) each edge, joining 2 sites on G, is replaced by a path of 3 edges. Therefore, 2 "additional" sites have to be introduced for each edge of G. The original sites are then "duplicated" along with their incident edges. For each of these pairs of of duplicated sites, one additional edge is added that connect the two sites of a pair. The resulting auxiliary graph G A is depicted in figure 3(b), where additional sites appear as squares and duplicated sites as circles. Figure 3(b) also illustrates the weight assignment on G A . A more extensive description can be found in [10].
(2) a MWPM on the auxiliary graph is determined [1]. A MWPM is a minimumweighted subset M of the edges contained in G A , such that each site of G A is met by one edge in M . This is illustrated in figure 3(c), where the bold edges represent M for the given weight assignment. The dashed edges are unmatched. Due to construction, the auxiliary graph consists of an even number of sites and since there are no isolated sites, it is guaranteed that a perfect matching exists.
(3) finally it is possible to find a relation between the matched edges M on G A and a configuration of negative-weighted loops C on G by tracing the steps of the transformation (1) back. Note that each edge contained in M that connects an additional site (square) to a duplicated site (circle) corresponds to an edge on G that is part of a loop, see figure 3(d). More precise, there are always two such edges in M that correspond to one edge on G. All the edges in M that connect like sites (i.e. duplicated-duplicated, or additional-additional) carry zero weight and do not contribute to a loop on G. Afterwards a depht-first search can be used to explore C and to determine the properties of the individual loops. For the weight assignment in figure 3(a), there is only one loop with weight ω L = −2 and length ℓ = 4. Note that the result of the calculation is a collection C of loops (and one path for (i), see below), such that the total loop weight is minimized. Hence, one obtains a global collective optimum of the system. Clearly, all loops that contribute to C possess a negative weight. Note that C can be empty and that sub paths are neither allowed to intersect nor to terminate at some site within the lattice.
So as to induce a MWP, as in problem (i), special care is needed. We have to allow the paths explicitly to terminate at a certain node. Therefore, two extra sites are introduced to the graph. One extra site is connected to each of the free boundaries by adding edges with weight zero to the transformed graph that join the site with the sites of the respecting boundary. Any subsequent MWPM will contain a path of minimal weight, joining the extra sites. This does not necessary coincide with the "shortest" path, as explained above, but yields the minimal-weighted path in the presence of negative weighted loops. In contrast to the loops, MWPs can carry a positive weight in principle, but, as we will see, this will not be the case close to and beyond the percolation transition we are interested in.

Minimal-weighted paths
Within this problem, the aforementioned collection C consists of a set of loops, which might be empty, and a MWP spanning the lattice between the free boundaries. Since the MWP is anyway of O(L) in the vertical direction by construction, we call the path percolating, if its projection on the horizontal axis (i.e. its roughness r) covers the full systems. We studied 2d lattices with ±J disorder and sampled over up to 3.5 × 10 4 (L = 128) realizations of the disorder to perform averages, subsequently denoted by . . . . Firstly, we investigate the percolation probability. For a very small fraction ρ of negative edge weights, the path will cross the lattice in a rather direct fashion, since overhangs are likely to increase the weight of the path and hence the weight of the whole configuration C. Increasing ρ also increases the spanning probability P x L (ρ), i.e. the probability of the horizontal extension to be O(L). It is expected to scale as P where the critical exponent ν describes the divergence of the correlation length at the critical point ρ c and f is a universal function. A data collapse, obtained using the method described in [9], yields ρ c and ν as listed in table 1 with a "quality" S = 0.79 of the scaling assumption. We want to point out that the value of ν found here is clearly distinct from the value 4/3 of conventional bond percolation, that would give a much worse quality of S = 5.20 (see also figure 4(b)). Note that many negative weights are needed to allow the MWP to percolate. And indeed, the path weight becomes negative above ρ ω c = 0.0869(2) with a probability that approaches unity quickly. Therein, the finite-size behavior is described by an exponent ν ω consistent with the value ν stated above, see figure 4.
Further note that P x L (ρ) < 1 even for large ρ, where the spanning probability seems to saturate slightly above 0.8. The reason is, as already mentioned, that we actually do not optimize the length of a single path but the weight of the whole configuration C. However, this behavior is clearly distinct from conventional percolation theory.
Next, we consider the probability P ∞ L ≡ ℓ /L d that a bond belongs to the path, where ℓ is the mean path length. It exhibits the finite-size scaling behavior where β signifies the percolation strength [13]. Here, we fixed ρ c and ν as obtained above and determined the value of β (S = 0.97) listed in table 1. Adjusting all three parameters in the above scaling assumption yields ρ c = 0.1027(1), ν = 1.45(4) and β = 1.06(2) with a quality S = 0.76. Note that the values agree within the errorbars.
We also determined some quantities (see also table 1) just at ρ c = 0.1032(5) with additional simulations up to L = 512 (2 × 10 3 samples). We studied the associated finite-size susceptibility χ L = L −d ( ℓ 2 − ℓ 2 ). Its finite size-scaling at ρ c can be described using the exponent γ via χ L ∼ L γ/ν . Furthermore, the mean path length shows the critical behavior ℓ ∼ L d f , where d f denotes the fractal exponent of the paths. Here, calculations at ρ c = 0.1027(1) would result in a slightly smaller value of d f . We further find the roughness exponent d r from r ∼ L dr to be compatible with unity. Note that the obtained exponents satisfy the scaling relations d f = d−β/ν and γ+2β = dν. We also measured the mean path weight ω p at the percolation point ρ c and found that ω p ∼ ℓ for L → ∞ seems to hold.
We can further probe C to not only investigate the MWP but to yield exponents that describe the small loops. A detailed study of the scaling behavior (length ℓ as function of the spanning length R, not shown here) at ρ c shows that they seem within error bars to exhibit the same fractal dimensiond f = d f as the MWPs. The loops also exhibit a mean (negative) weight that increases linearly with the loop length ℓ. Further, one expects the distribution n ℓ of loop lengths ℓ to exhibit an algebraic scaling n ℓ ∼ ℓ −τ , where τ signifies the Fisher exponent. Here, finite-size effects and the presence of the path lead to a suppression of large loops and thus to a deviation from the expected scaling behavior already for rather small values of ℓ. In this question, we find a most reliable value τ if we account for corrections to scaling via n ℓ ∼ ℓ −τ /(1 + bℓ ω ), see table 1 and inset of figure 6. In principle, at criticality, the Fisher exponent and fractal dimension are related through the scaling relation τ −1 = d/d f . This would lead us to expect τ ≈ 2.58.
Finally, in 2d we can associate a cluster with each loop. The respective volume v is measured as the number of enclosed plaquettes. It scales as v ∼ Rd v with d v = 2.00(1), revealing the compact nature of the loops interior.

Minimal-weighted loop configurations
Here, we studied 2d (hex, 3d) lattices with size up to L = 128 (192, 48) and 4.5×10 4 (1.2×10 4 , 1.6×10 4 ) samples. First, we analyze the largest loop for a given realization of the disorder. We find the linear extensions of the loops by projecting it onto all perpendicular axes. A loop is said to span the lattice, if its projection completely covers at least one axis. From the probability P s L (ρ) that a loop spans the lattice we estimate ρ c and ν listed in table 1. The qualities of the scaling law were found to be S ± 2d = 1.09, S ± hex = 0.79, S Gl 2d = 0.91 and S ± 3d = 1.9. Further, it is interesting to note that the mean number of spanning loops N satisfies a similar scaling relation as the percolation probability, see figure 5(a), governed by the same values for ρ c and ν as P s L (ρ). In either case we find N > 1 for large values of ρ, as mentioned already in the introduction, see also [11]. As above, the probability P ∞ L (ρ) that an edge belongs to the percolating loop can be used to determine the exponent β, see figure 5(b) and table 1. Herein, the qualities of the scaling law were S ± 2d = 1.88, S ± hex = 0.32, S Gl 2d = 1.16 and S ± 3d = 2.08. Interestingly, the values of ρ c , ν and β for the 2d loops with ±J disorder are reasonable close to those of the 2d random-bond Ising model with ρ c = 0.103(1), ν = 1.55(1) and β = 0.9(1), see [3] and references therein. This probably means that the 2d ferromagnet to spin glass transition at T = 0 can be explained in terms of a percolation transition in the following way: For a spin glass, one starts with a ferromagnetic configuration and searches for loops in the dual lattice with negative weight. These loops correspond to clusters of spins, which can be flipped to decrease the energy. If these loops are small, i.e. not percolating, the ground state is ferromagnetic, otherwise the ground state exhibits spin-glass order. A different argument relating this transition to a percolation transition was also briefly mentioned in a study, which focuses on the critical slowing down of polynomial-time algorithms at T = 0 phase transitions [17].
Right at ρ c we studied 2d (hex, 3d) systems up to L = 512 (768, 96) with 3.2×10 4 (1.6 × 10 3 , 6.4 × 10 3 ) samples. For the 2d systems we found γ listed in table 1 with  Figure 6. Results for the 2d ±J loops: distribution n ℓ of the loop-lengths ℓ at criticality, excluding the spanning loops (main plot), gray data points were omitted from the fit. The inset shows n ℓ for the case of MWPs, where the fit accounts for corrections to scaling (see text).
qualities Q ± 2d = 0.61, Q ± hex = 0.10 and Q Gl 2d = 0.51. In 3d, due to the very small percolation probability 0.0014(1) at ρ c our results are less clear. Hence, we can draw no final conclusions for that case. Regarding the fractal dimension, spanning loops and the non-spanning loops exhibit within error bars the same exponents d f =d f in 2d (values of d f listed in table 1). As in the MWP case, the non-spanning loops are compact ( v ∼ R 2 ). On the other hand, in 3d, we findd f = 1.43(2) different from d f = 1.30 (1).
Also regarding the mean weight ω ℓ for given length ℓ, we find in 2d that spanning and non-spanning loops exhibit both ω ℓ ∼ ℓ (i.e. ∼ L d f for the spanning loops), while in 3d the quality of the data is again not sufficient to observe a clear power law.
The distribution of loop-lengths, excluding the truly spanning loops, is in 2d in good agreement with a power-law decay n ℓ ∼ ℓ −τ , see figure 6 and table 1. In 3d we find again strong finite size effects and hence a most reliable estimate of τ using a scaling form n ℓ ∼ ℓ −τ /(1 + bℓ ω ) .
Finally, we address the performance of the algorithm. We doe not want to measure the running time in terms of CPU minutes, since this is machine dependent and is also influenced by external factors like which other processes are running. Hence, we have to look at the algorithm more closely. So as to determine a minimumweight perfect matching, the algorithm attempts to find an optimal solution to an associated dual problem [6,8]. Therein, the basic task of the algorithm is to find augmenting paths that improve the solution to the latter. While executing, the solution to the dual problem is adjusted several times. As a measure for the algorithm performance we consider the number N da of such adjustment operations, until the algorithm terminates. Figure 7 shows the average number of these operations per lattice site n da = L −d N da for the case of a Gaussian-like distribution of the edge weights. The value of n da increases with increasing ρ and curves for different system sizes deviate from each other as the critical point is approached from below. At ρ c it scales like n da ∼ L 2.074 (4) . Moreover, the corresponding susceptibility χ da , i.e. the fluctuations of the n da , diverges at the critical point where is scales as χ da ∼ L 1.41(4) , as can be seen from the inset of figure 7. For the case of a ±J distribution of the edge weights, we found similar but slightly different results (not depicted): Here, data sets that describe the scaling of n da for different system sizes, fall onto one universal curve everywhere. This holds also for χ da , which is peaked slightly below the critical point at ρ ≈ 0.096.

Conclusions
We have introduced a percolation problem, where edge weights with possibly negative values are attached to the edges. A system is called percolating if spanning loops or paths of negative weight exist. Hence, this model might be fundamentally different from classical percolation problems. We studied systems in 2d and 3d up to large system sizes. In all cases, the universality class is indeed clearly different from standard bond/site percolation. Further, the values of the fractal exponents obtained for the paths/loops differ from those reported in the context of shortest paths or hulls of percolation clusters, so far. Nevertheless, the numerical values for the exponents found here are close to those of disorder-induced single defects (d f = 1.261 (16)) and multiple defects (d f = 1.250(3)) in a 2d elastic medium at zero temperature [16].
We observe universality in 2d, hence the properties are independent of the type of weight distribution, lattice geometry and the same for loops or paths. We also studied (not shown here) a related problem, where the loops are allowed to intersect. The corresponding mapping was recently used in a different context to determine extended ground states of Ising spin glasses [15]. We found again the same critical exponents. On the other hand, we observe in 2d the same behavior as for the T = 0 ferromagnet to spin glass transition for the random-bond Ising model, hence this physical transition can be probably explained by a percolation transition in the spirit of the model we have introduced here. Hence, studying percolation problems with negative weights and similar generalizations might be a key approach to describe many yet not wellunderstood phase transitions in terms of percolation transitions. We also anticipate applications to other fields like social problems (see introduction) or other types of networks.