Biased and greedy random walks on twodimensional lattices with quenched randomness: the “greedy” ant within a disordered environment
Abstract
The principle characteristics of biased greedy random walks (BGRWs) on twodimensional lattices with realvalued quenched disorder on the lattice edges are studied. Here, the disorder allows for negative edgeweights. In previous studies, considering the negativeweight percolation (NWP) problem, this was shown to change the universality class of the existing, static percolation transition. In the presented study, four different types of BGRWs and an algorithm based on the ant colony optimization (ACO) heuristic were considered. Regarding the BGRWs, the precise configurations of the lattice walks constructed during the numerical simulations were influenced by two parameters: a disorder parameter that controls the amount of negative edge weights on the lattice and a bias strength that governs the drift of the walkers along a certain lattice direction. Here, the pivotal observable is the probability that, after termination, a lattice walk exhibits a total negative weight, which is here considered as percolating. The behavior of this observable as function of for different bias strengths is put under scrutiny. Upon tuning , the probability to find such a feasible lattice walk increases from zero to one. This is the key feature of the percolation transition in the NWP model. Here, we address the question how well the transition point , resulting from numerically exact and “static” simulations in terms of the NWP model can be resolved using simple dynamic algorithms that have only local information available, one of the basic questions in the physics of glassy systems.
pacs:
05.10.Ln, 05.40.Fb, 05.70.JkI Introduction
The purely stochastic motion displayed by a simple random walk (RW) and its various extensions and modifications that, e.g., go under the popular synonyms “the ant in a labyrinth” de Gennes (1971); Straley (1980) and “the biased ant in a labyrinth” Stauffer (1999), give rise to intriguing effects such as (anomalous) diffusion, drift and trapping. These, and further related topics, were actively studied since the beginning of the last century, first by analytic means and later using extensive numerical simulations. A subset of those studies, relevant in the context of the presented article, is reviewed in subsect. II. In all of these studies, a random walker is considered that traverses a given (disordered) environment subject to different dynamicsgoverning rules.
Here, we also consider a random walker on a disordered “energy” landscape, where, in contrast to most previous studies, the underlying energyvalues might have a positive or a negative sign. Note that the existence of the negative energies fundamentally changes the behavior of the model defined below and will be destroyed by lifting all energies by the same amount to a positive value. The reason for this fundamental difference is that, in the presence of a suitable amount of (nonreplenishable) negative energies, the walker might lower its energy more and more by performing longer and longer walks, in contrast to the presence of only positive energies. The energy values are drawn from specified distributions wherein a disorderparameter allows to alter the fraction of negative energies in the environment (see discussion below).
Regarding such a disordered environment, we have previously introduced Melchert and Hartmann (2008) and investigated Apolo et al. (2009); Melchert et al. (2010, 2011); Norrenbrock et al. (2012); Claussen et al. (2012) the negativeweight percolation (NWP) problem by numeric means. Regardless of the spacial dimension of the underlying (hypercubic) lattice graph, the observables in the NWP problem are always linelike, i.e. have an intrinsic dimension of . As observables in the NWP problem one considers, e.g., minimumweight paths (in the presence of a possibly empty set of negative weighted loops) that span the underlying lattice between a pair of specified boundaries. The problem of finding these paths numerically can be cast into a minimumweight perfect matching problem, e.g. discussed in Refs. Melchert et al. (2010); Norrenbrock et al. (2012) in more detail (in the latter reference, a comprehensive description and illustration of the respective algorithm can be found). A pivotal observation is that, as a function of the disorder parameter , the NWP model features a disorderdriven phase transition, at which the actual minimumweight path energy becomes negative. In this question, note that the aforementioned transition in the path energy is a precursor to a further “geometric” phase transition at which the paths exhibit a roughness of the order of the systems size and where loops might appear that span the system along at least one direction Melchert and Hartmann (2008). In the limit of large system sizes, there is a particular value of the disorder parameter, signified as , at which paths with a negative energy appear for the first time. E.g. for the square () lattice, considering a bimodal disorder distribution, the average path energy turns negative above Melchert and Hartmann (2008). Therein, the finitesize scaling behavior is described by the critical exponent , consistent with the one that describes the geometric transition of the loopsonly setup. At , slightly above this critical threshold, the fractal dimension of the paths was estimated to be . Further, at this critical point, the loops that one finds in addition to the minimumweight path are rather small and isolated. Hence, they (presumably) do not affect the precise configuration of the minimumweight paths.
As mentioned above, NWP was previously studied numerically exact by employing a mapping to an auxiliary minimumweighted perfect matching problem, which relies on an involved optimization algorithm specifically tailored for the problem at hand. Instead of using such a complex algorithmic approach that finds the optimum of a given optimization problem using global information, we here address the question how well the aforementioned transition point (resulting from numerically exact simulations that involve a high degree of optimization) can be resolved using simple dynamic algorithms that have only local information about the assignment of edge weights available. This is a more physical viewpoint on the problem, where, typically, the movement of particles is mostly determined by local information. Hence, a comparison between global static and dynamical physical viewpoints is possible. These two perspectives are often taken for the study of disordered and glassy systems Young (1998); Hartmann and Weigt (2005); Binder and Kob (2011). Therefore we consider different biased random walk approaches and study their ability to find system spanning paths with negative pathweight. Averaging over different realizations of the disorder, we compute the probability to find a negative weighted path as a function of the model parameter for different simple dynamics. The algorithms we consider here all mimic moving particles (also referred to as agents or walkers) in a disordered environment. The simplest dynamics we study stems from a biased random walk on the lattice and the most intricate dynamics is implemented in terms of a particular heuristic algorithm known as “ant colony optimization”. Subsequently, we will refer to as the “dynamic” transition points, as opposed to the result from the “static” (global) simulations quoted above. Not too earnest one might allude to this study as being dedicated to the behavior of “the greedy ant in a disordered environment”.
The remainder of the presented article is organized as follows. In section II we review some of the related literature that alludes to effects that we also expect to observe within our simulations. In section III, we introduce the model in more detail and we outline the different dynamic algorithms used to perform the lattice walks. In section IV, we list the results of our numerical simulations and in section V we conclude with a summary.
Ii Review of related literature
A typical observable for the characterization of ordinary RWs is given by their mean square displacement (MSD) , which, as a function of time , is expected to scale as . Therein, signifies the diffusion coefficient characteristic for the observed walks.
A particular extension of the simple RW model, relevant in the context of the presented study, consists in the addition of disorder to the environment in which a walk is performed. Different approaches to implement disorder can be found in the literature: E.g., one might consider percolationlike “occupied vs. free” site (or bond) disorder BenAvraham and Havlin (1982); Gefen et al. (1983); Pandey (1984a); Argyrakis and Kopelman (1984); White and Barma (1984); Havlin and BenAvraham (1987), in which case a fraction of lattice sites (or bonds) is distinguished as “occupied”, and where a RW is performed on a cluster of occupied nearest neighbor sites (a.k.a. “the ant in the labyrinth”). For short walking times, i.e. times within which the walker cannot fully trace the cluster, and in the vicinity of the critical point , signaling the onset of percolation, anomalous diffusion was observed Gefen et al. (1983); Bouchaud and Georges (1990). Thereby, the MSD of the RW scales according to , characterized by an effective dimension . At the critical point for site dilution on lattices, Ref. BenAvraham and Havlin (1982) reports , i.e. indicating a subdiffusive behavior. For the non disordered case (i.e. in the limit ) one would expect to find , characterizing the effective “fractal” scaling dimension of unhindered RWs.
A different approach considers “energetic” disorder associated with sites or bonds Argyrakis et al. (1984); Avramov et al. (1993), representing random barriers that define sitetosite transition probabilities governed by Boltzmann statistics. Therein, a model inherent temperatureparameter controls the ability of an individual walker to overcome energy barriers. At low temperatures, a RW is likely to get trapped at local minima in the energy landscape. Albeit this trapping effect is only temporary, the RWs generally display a rather limited mobility. With increasing temperature the RWs gain mobility and behave similar to ordinary RWs in the limit . Regarding this energetic disorder one might further distinguish static (i.e. quenched) disorder, which remains unchanged in time, and dynamic disorder, which is renewed after each discrete time step, see Ref. Avramov et al. (1993). In case of static disorder and after a temperature dependent crossover timescale , the MSD exhibits the same asymptotic linear scaling as for unhindered RWs. Supported by analytical arguments from classical chainreaction theory, the crossover time could be attributed to an effective activation energy that must be overcome in order to escape from an initial local minimum in the energy landscape. Intuitively, the lower the temperature, the larger the crossover time tends to be. In the early time regime, i.e. at times smaller than , the scaling of the mean square displacement is subdiffusive. In case of dynamic disorder, no crossover time is found and the same scaling as for the unhindered RW is observed for all temperatures. The renewal of the disorder after each time step creates an averaged environment in which the temperature directly reflects the transport efficiency in the medium, resulting only in a temperature dependent diffusion constant .
A further extension that is heavily studied in the literature consists in emphasizing a particular lattice direction along which the RW effectively drifts. Also for this modification there are different possible implementations, differing in the particular way in which the transition probabilities are altered, see Bunde et al. (1987); Arapaki et al. (1997); Avramov et al. (1998); Dhar and Stauffer (1998); Stauffer (1999); Kirsch (1999); Zhang et al. (2006). In this regard we only briefly recap one of those studies, which, for the first time, reported on a direct observation of a sharp drift to nodrift transition Dhar and Stauffer (1998). Therein the influence of a directional bias field on the diffusive behavior of random walks on lattice graphs with percolationlike site disorder at , i.e. above the critical point of the setup, was addressed. Initially, for a given realization of the disorder, noninteracting RWs were started on randomly chosen occupied sites and individual walkers were only allowed to advance to adjacent occupied sites. A bias parameter is used to control the drift of the RWs along the xdirection: the probability to advance along the positive xdirection is set to , the probability to advance to any of the six neighbors is set to . Once such a step is proposed, it is only executed if the target site also belongs to the cluster. Hence, there is the possibility that a walker remains at its current position (speaking in terms of the popular phrases mentioned in the introduction, it would be more precise to call this the “blind biased ant in the labyrinth” since the decision on the direction of the proposed move is made irrespective of whether the targetsite can actually be reached). For such noninteracting walkers and prior to the numerical study reported in Dhar and Stauffer (1998), a drift to nodrift transition at a finite critical value was theoretically hypothesized M. and Dhar (1983); Dhar (1984). For weak bias , the walkers exhibit a diffusive motion and their velocity assumes a constant value which tends to increase with . For strong bias, diffusion is slowed down since the walkers tend to wind up in “dangling ends” of the clusters on which they reside. With increasing bias they have an increasingly hard time to escape from such dangling ends. Neglecting initial transients, three different regimes were observed: (i) , where tends to a constant, (ii) , where , and, (iii) , where (and ). The critical bias strength at which the drift to nodrift transition occurs was estimated as . The behavior of the average velocity for a large range of values, covering all three scaling regimes, could further be summarized by an extrapolated scaling form. Similar studies were, e.g., also carried out to construct a full phase diagram for the drift to nodrift transition in the –plane for the square lattice Kirsch (1999). Regarding the MSD, the RW component of the dynamics leads to diffusion and the bias induces a drift . At this point note that under drift, the asymptotic effective dimension might be expected. At short (long) times diffusion (drift) is dominant. Under a strong bias, the dangling ends act as traps that inhibit drift by characteristic “trapping times” (depending on the precise cluster structure). Consequently, the MSD is governed by effective exponents that approach their asymptotic value rather slowly Pandey (1984b); Dhar (1984); Kirsch (1999); Stauffer (1999). Further, biased diffusion for networks without underlying regular spatial structure (i.e. regular random graphs, ErdősRenyi random graphs and scalefree networks) was recently investigated by means of numerical simulations and analytic calculations Skarpalezos et al. (2013). Therein, the bias was designed in a way that it guides the walker to a specified target node along the respective shortest path. Using a model parameter, the strength of the bias, ranging from unbiased (yielding an ordinary random walk dynamics) to fully biased (effectively confining the walker to the shortest path), could be tuned. Among other things it was shown that for regular random graphs and ErdősRenyi random graphs, the scaling of the mean first passage time as function of the system size changes from a powerlaw to a logarithmic behavior at a characteristic value of the bias parameter.
In addition to simple RWs, selfavoiding walks (SAWs), conveniently used as models for randomly bent polymers that exhibit an “excluded volume” effect Caracciolo et al. (1992), have been studied for diluted systems exhibiting percolation disorder Kremer (1981); Grassberger (1993). Therein, a basic observable is the mean square radius of step SAWs which scales according to , where the exponent can be seen as an inverse of a SAW characteristic scaling exponent that is universal in each dimension . Depending on the fraction of occupied sites on the underlying lattice and the set of configurations over which averages are computed, different SAW scaling regimes can be found Grassberger (1993). For the case of no disorder (equivalent to ), results consistent with () are obtained in two dimensions, see, e.g., Ref. Berretti and Sokal (1985). A different study, aimed at clarifying the geometric and energetic scaling behavior of minimum energy (ME) SAWs on lattices with random site energies Smailer et al. (1993), reports () for the case (obtained via exact enumeration methods). Thus, under the influence of quenched randomness, minimum energy SAWs tend to expand. Note that within errorbars this value is consistent with the fractal dimension of paths in the NWP problem at the particular point (also in the estimates of the respective exponents stack up well: as compared to Melchert et al. (2010)).
A further variant of a simple RW, wherein a loop is erased as soon as it is formed, is referred to as looperased random walk (LERW) Majumdar (1992). A LERW might be interpreted as a simplified version of a SAW which has some correspondence to spanning trees and, in , exhibits a scaling dimension of Grassberger (2009); Majumdar (1992).
Iii Model and algorithms
In the present article we consider lattice graphs with a square lattice geometry, having side length and periodic (free) boundary conditions (BCs) in the vertical (horizontal) direction. The considered graphs have nodes (plus two extra “outer” nodes) and undirected edges that join adjacent nodes on the lattice, see Fig. 1. In addition, there are extra edges that attach the nodes on the left and right boundary to the respective outer nodes. We further assign a weight to each . These weights represent quenched random variables that introduce disorder to the lattice. Here we consider independent identically distributed weights drawn from (i) the bimodal distribution
(1) 
where signifies the fraction of negative edge weights , or (ii) the semicontinuous distribution
(2) 
wherein signifies a random number uniformly drawn from the interval .
Subsequently we will study lattice walks on instances of weighted lattice graphs with sidelength up to . The considered walks are confined to start at an outer node attached to the left lattice boundary and end as soon as the walker reaches an outer node attached to the right boundary, see Fig. 1. More precisely, to construct the lattice walks, two classes of algorithms were used: biased and greedy RWs and an ant colony optimization (ACO) heuristic.
iii.1 Computing paths via biased and greedy random walks (BGRWs)
In brief, the biased RWs explore the lattice graph while greedily choosing the “best” direction at each discrete timestep. This might be the lowest weighted edge connected to the node the walker currently resides on, or, if there is a draw, one of the lowest weighted edges picked at random. By “biased RW” it is meant that following a certain rule, a preferred direction is chosen for the next step, even if it is not the best possible immediate choice. The rule employed during our simulations reads as follows: say the preferred direction is given by the positive xdirection. First, the walker determines the direction corresponding to the best local edgechoice. If this direction is already the positive xdirection, the walker goes there. Otherwise, if the best local choice does not lead towards the positive xdirection, the walker accepts the choice only with probability , where (termed “bias probability” or short “bias”). I.e., if its not the best local choice, the walker goes towards the positive xdirection with probability .
Thus, a walker experiences a drift that, depending on the biasdirection, effectively “guides” him towards a particular direction (see discussion in sect. II). Here, we will refer to this kind of RW in a disordered environment as “biased and greedy random walk” (BGRW). While traversing the graph, the walker keeps track of the edge weights it encounters. Consequently, a lattice walk is considered successful if the total sum of the edges traversed by the walker is smaller or equal to zero. Clearly, the walker will be able to find negativeweight paths only if they exist, which can be detected by the previously mentioned exact algorithms. We anticipate that there will be an intermediate range, where negativeweight paths exist, but are so rare that they cannot be detected by the locallyacting walker.
We also study a variant of such lattice walks, where the walker interacts with the graph by modifying the weights of the traversed edges. This is to some extent motivated by studies, as, e.g., the one reported in Ref. Boyer et al. (2004), aimed at mimicking the foraging behavior of social monkeys. Here, we chose the rules that the walker

replaces both, positive and negative edge weights by weight (aimed at modeling a finite nonreplenishable amount of both, resources and costs), or

only replaces negative edge weights with the standard positive value of (aimed at modeling nonreplenishable finite resources, only).
These types of modifications where chosen to also prevent the walker from getting stuck in regions of the lattice with a high density of negative edge weights, as for RWs in disordered energy landscapes at low temperatures reviewed in sect. II above. If no such interaction with the lattice is implemented, the walker could in principle collect unlimited amounts of negative edge weights.
Four distinct variants of BGRWs, termed RWs of type A, B, C or D, were considered to evaluate which one performs better in finding paths with an overall negative weight. Here, the term “better” means “with higher probability for a given disorder parameter value ”. Also bear in mind that ultimately, the aim is to quantify to which extent the transition point (known from numerical simulations that involve finding the exact global optimum Melchert and Hartmann (2008)) can be resolved using simple dynamic algorithms that have only local information about the assignment of edge weights available. The four variant of biased RWs, along with a brief description of the dynamic governing rules, are listed below. The variant A is based on standard RWs, while the other three variants are based on (partially) loop erased RWs.

BGRWs of type A: this type of walk is given by a BGRW, that replaces negative edge weights (rule (ii) above) and maintains a running sum of the edge weights encountered while traversing the graph.
For completeness, note that also the variant wherein the walker replaces edge weights according to rule (i) was considered, see Fig. 2. In the figure, type A walks on a square lattice of size for different values of the disorder parameter are shown, considering the bimodal disorder distribution (see Eq. 1) for different fractions of negative edge weights .
While traversing the lattice graph, the walker modifies positive and negative edge weights, replacing them by edge weights with value zero. As a consequence, for a small fraction of negative edge weights (as for in Fig. 2(a)) it can be observed that the walker backtracks along subpaths of considerable length. This is due to the fact that, after the walker modified the respective edge weights, they are actually “cheaper” than the edge weights that characterize the majority of the edges. The bias is designed to induce a drift along the positive xdirection, hence one finds a pronounced “channeling effect” along the horizontal lattice axis. Once such a region of zero weighted edges develops, for small values of , the walker has a hard time to escape from it. As evident from Figs. 2(b,c), the channeling effect is less pronounced as the fraction of negative edge weights increases. Note that this effect is conceptually similar to the behavior of RWs in energetically disordered lattices, reviewed in sect. II (see Refs. Argyrakis et al. (1984); Avramov et al. (1993)), where initially an effective activation energy must be overcome in order to escape from a local minimum in the energy landscape. The general characteristics of such a lattice walk, basically confined to take place on edges of zero weight in an environment that mostly features positive edge weights, further bears some resemblance to the trapping behavior observed for RWs on percolation clusters discussed previously.

BGRWs of type B: this type of walk is given by a looperased BGRW, traversing the graph without modifying the edge weights (which is not necessary, due to the erasure of the loops). The sum of edge weights along the looperased walk is computed as soon as the BGRW has finished.

BGRWs of type C: this type of walk combines an interaction with the environment (as for type A walks) with looperasure (as for type B walks). The lattice walk is constructed via the following twostep procedure: In the first step, the walker performs a BGRW and modifies the edge weights related to the traversed edges according to rule (i) above. At each visited node it further stores the sum of edge weights accumulated so far, the previous node it has visited, and the original weight (i.e. the value of the edge weight before it was modified by the walker) of the edge it has just traversed. In the second step, after the BGRW has terminated, the looperasure and subsequent path evaluation is performed. Therefore, the lattice path is traced back, starting from the final node, by selecting at each intermediate node the step that yields the lowest sum of the edge weights until that point. This selection appears each time when a loop is encountered while backtracking. This is a heuristic which leads to lower path weights, but is not able take combinatorial cases like a loops inside loops inside loops, etc., into account, since this would contrast the idea of using fast and local heuristic algorithms for computing the walks. Finally, the total weight of the path is computed by summing up the edges. For completeness, note that also the variant wherein the walker replaces edge weights according to rule (i) was considered.

BGRWs of type D: this type of walk might be referred to as “partially” loop erased BGRW, wherein only unfavorable loops along the lattice walk are erased. More precisely, the walker modifies the weights of the traversed edges and checks each time if the node it currently resides on was visited previously. If so, the walker checks whether the created loop has a negative or positive weight. The loop is kept as part of the lattice walk if it has a negative weight, or, if it has a positive weight, the loop is discarded and the original edge weights are restored along it.
Regarding these dynamics, only the variant where the walker replaces edge weights according to rule (ii), i.e., the modification of negative edge weights only, was considered. This was done since the modification of positive and negative weights led to the walker being stuck for a very long time in a certain part of the graph, especially when using a low value for the bias and a small fraction of negative weights. The cause for this was that the walker created areas of zero weight in the graph that it could not escape efficiently if the surrounding edges were all positive and the bias was too low (somehow similar to the effective activation energy that must be overcome in the initial phase of RWs in energetically disordered lattices reviewed in sect. II). A similar horizontal “channeling” effect can be seen in Fig. 2 for random walks of type A BGRWs, where the walker is constrained to move on an almost straight line because of the “repulsive” effect of the positive edge weights.
iii.2 Computing paths via ant colony optimization (ACO)
Previously we considered algorithms that took into account local information only. I.e. the random walk approaches we considered so far made a decision on where to go next by only considering the weights of the edges adjacent to the walkers current location.
In this subsection we will consider a nature inspired, population based heuristic for the solution of optimization problems, called Ant Colony Optimization (ACO) ACO (); Dorigo and Blum (2005). Designed to mimic the foraging behavior of actual ants, the ACO heuristic is a particular example of swarm intelligent systems, which has proven to be valuable in solving a wide range of optimization problems that can be cast into the form of optimalpath problems. As regards this, the ACO heuristic has already been applied to different optimization problems such as, e.g., the traveling salesperson problem Dorigo (1997), which is a hard optimization problem, and the minimumweight spanning tree problem Neumann and Witt (2010), which is polynomially solvable, i.e., easy.
The ant colony represents a population of, say, individual agents that are all able to construct solutions to the given problem by considering local information only (these need not necessarily be solutions of high quality). Albeit the agents do not interact directly, they are able to interact indirectly through the deposition of pheromone on the graph edges (see discussion below). These interactions should lead to lowerenergy solutions as compared to the single randomwalk algorithms discussed above.
Here, the problem is to find a loopless minimum weight path from a specified source node to a specified target node in a weighted undirected graph , as discussed earlier. The ACO heuristic comprises an iterative algorithm, where in the beginning, there is no pheromone on the edges . Hence, the only information associated with a particular edge is its weight . A basic variant of the ACO algorithm for the above problem can be cast into the following steps:

Preprocessing of the edge weights: Transform the set of edgeweights to a set of (initial) transition probabilities , where with specifying the set of nodes adjacent to node . In this way it is assured that an edge weight with a comparatively large negative weight will result in a comparatively high probability for the agent to take the step . These initial transition probabilities can also be thought of as an initial distribution of pheromone which affect the behavior of the individual agents.

Generate solution to the problem: The ability of an individual agent to construct loopless paths as solution to the problem at hand can be implemented in various ways. For a more clear illustration lets consider one agent only, i.e. . Here, starting at a specified source node , we let the agent perform a loop erased random walk (LERW), guided by the transition probabilities . During the loop erased walk, the agent sums up the edge weights along the LERW edges. As soon as the agent reaches a specified target node , the guided LERW is completed and a solution to the problem, i.e. a loopless path, is obtained. Note that this is not necessarily a solution with a good quality regarding the given optimization criterion.

Evaluate quality of the solution and deposit pheromone: In order to quantify the quality of the LERW obtained in step (ii) and so as to reflect the optimization criterion of the problem at hand, we compute the “fitness” parameter , where is a tunable parameter, while and are the weight and length of the LERW, respectively. Note that for edge weights drawn from a bimodal distribution, it holds that . Further, note that the more negative the sum of the edge weights, the larger the value of is. Then, modify the transition probabilities along the LERW to and normalize them as in step (i) above. This is called a delayed pheromone update, since the local information on the edges is updated after the LERW is obtained. The parameter can be tuned to alter the influence of the quality measure on the subsequent recomputation of the transition probabilities. Here, we optimize the ratio , effectively measuring the average edgeweight that contributes to the path. This quantity was chosen since it naturally normalizes the value of to the range and thus allow to control the influence of the current local solution on the transition probabilities in an easy way by using the parameter (in contrast, note that in the static NWP problem we optimize instead).

Evaporate pheromone: While pheromone is deposited only on the edges of particular loopless paths, it can be useful to also evaporate some of the pheromone from all edges. E.g., a simple evaporation heuristic is to modify the transition probabilities along all edges to , where , and to normalize them as in step (i) above.
Step (i) is a preprocessing step to generate transition probabilities that reflect the weight assignment on . Steps (ii) and (iii) complete the “lifecycle” of an individual agent: It generates a (nonoptimal) solution to the optimization problem and evaluates a fitness for the solution. The fitness is then used to alter the transition probabilities which are associated with the edges. Instead of using just one agent, steps (ii) and (iii) straightforwardly generalize to individual agents. The construction and evaluation of LERWs and the subsequent pheromone evaporation, i.e. steps (iiiv) comprise one sweep. Consequently, the ACO algorithm is iterated for a number of, say, sweeps. The transition probabilities change from sweep to sweep. Thereby, edges that belong to paths with a high fitness get equipped with a larger transition probability. Albeit this induces a positive feedback that allows to distinguish “good” paths, the parameter can be used to limit the extent to which the transition probabilities are modified from sweep to sweep (during our simulation we used ). Also note that the parameter serves to limit the “longterm memory” of the search process in order to efficiently explore a large variety of paths (during most of our simulations we used ). E.g., setting and completely suppresses the deposition and evaporation of pheromone. For that parameter choice the individual agents would construct LERWs which are guided only by the transition probabilities induced by the edgeweights on the lattice (following step (i)) and they would not be able to indirectly communicate through the deposition of pheromone. Also note that in the extreme case where and , the edges that belong to a path found in step (ii) would get a transition probability and all transition probabilities that relate to those edges that depart from that particular path would effectively be zero. Hence, once an agent crosses that path (during a later sweep), it is immediately “confined” to that path and is unlikely to escape from it. During a simulation run for a particular realization of the edge weights we store the best path found so far. This path is returned by the ACO algorithm after it terminates.
Iv Results
As pointed out in the introduction, minimum weight paths on disordered lattices were already studied in terms of the NWP model, which, for a given realization of the disorder, can be solved numerically exact by using quite involved optimization algorithms. Instead of using such a complex algorithmic approach that takes into account global information we here address the question how well local algorithms perform. In particular we investigate how well the transition point for the square lattice with bimodal disorder, resulting from static simulations considering the NWP model, can be resolved using simple dynamic algorithms. We study local algorithms because the movement of particles in physical systems, or the decisions of agents in populations, typically depend on local information only. Hence, such an approach reflects real situations better. Therefore we consider the dynamic biased and greedy RW dynamics discussed previously in Sect. III.1 and the ant colony optimization heuristic described in Sect. III.2. The results are reported in Sects. IV.1 and IV.2, respectively.
iv.1 Biased and greedy random walks
At first, we consider type A BGRWs where the walker interacts with the environment via rule (i), i.e. while traversing the graph it replaces positive and negative edge weights by edge weights equal to zero. In Fig. 3 we show the probability that, after termination, a lattice walk with negative weight is found. Simulations are carried out for systems of size where data points represent an average over realizations of the disorder.
To facilitate intuition, for a bias strength the greediness of the walker has no effect. It starts at the outer node on the, say, left lattice boundary and performs a straight line (i.e. step) walk to arrive at the outer node on the right lattice boundary (thus, the lattice walk is characterized by a fractal dimension ). Since the disorder is uncorrelated, the walker basically sums up uncorrelated random edge weights drawn from the disorder distribution , see Eq. (1). At the value of the disorder parameter and for large values of (on average) half of the edge weights will be negative. Hence, above , the value of will tend to as . At this is equivalent to a symmetric RW, where, starting at the origin, one allows the walker to perform steps, afterwards asking for the probability that the walker is located, say, left of the origin. Also note that in this limit, the effective problem statement is trivially equivalent to the NWP problem. To further understand the scaling behavior of the data curves for the case shown in Fig. 3, one might follow a realspace renormalization approach Wilson (1971, 1983); Stanley (1999). The basic idea is to replace a subsystem of edges (characterized by disorder parameter , which is, generally spoken, the probability to have a negative weight) by one edge characterized via an effective disorder parameter . In this regard, note that the only condition to yield a “successful” negative weighted step path is that the number of negative edge weights along the lattice walk needs to exceed the number of positive edge weights by at least one. In this question, the probability to find a negative weighted path for systems of odd for a given value of is simply a sum of binomial terms:
(3) 
where . The fixed points for , which charazterize phases/ phase transitions in the renormalization, are determined by . Intuitively, a stability analysis of the three fixed points reveals that the only unstable fixed point is the previously discussed value (at ). Linearizing about this critical point yields the scaling power of the observable Eq. (3) via
(4) 
Albeit the expression for converges rather slowly, e.g. it assumes the values , and , it attains the limiting value as (obtained through a series expansion of the derivative of Eq. (3) in powers of ). This slow convergence is presumably caused by the additive character of edge weights: the total sum of the edges must be negative, while subsystems of arbitrary size are free to have a positive weight. Consequently, if the data curves for obtained in the simulations exhibit the property of scaling it should thus not come as a surprise if the scaling behavior is governed by the scaling parameter . In deed, we find that the data curves show a smooth data collapse if rescaled according to
(5) 
where for the parameter values and yield a good data collapse (not shown). This scaling behavior still holds for values of , where . E.g., in Fig. 3 the scaled data for is illustrated, where and .
BGRW type  (bi)  (bi)  (sc)  (sc) 

A (rule (ii))  0.30(1)  0.344(4)  0.27(1)  0.464(5) 
A (rule (i))  0.18(1)  0.199(4)  0.12(1)  0.280(3) 
B  0.0  0.258  0.0  0.465 
C (rule (ii))  0.23(1)  0.300(4)  0.20(1)  0.401(6) 
C (rule (i))  0.18(1)  0.174(6)  0.12(1)  0.241(4) 
D  0.0  0.187  0.0  0.286 
For decreasing bias strength the greediness of the BGRW has an increasing influence on the dynamics of the walker. Hence, at a given value of the average path weight after a fixed number of steps is smaller at smaller bias strengths. Consequently the probability to find a negative weighted path principally increases. Albeit true for intermediate bias strengths , this does not hold for comparatively small values . As evident from Fig. 4, the critical point above which negative paths appear first as tends to increase for a decreasing bias strength . I.e., one can observe a minimum value of at a characteristic value (for type A BGRWs considering an interaction with the environment according to rule (i)), see Tab. 1. A further extremal parameter set for type A BGRWs is the case of vanishing bias strength at . Since indicates the onset of a system spanning (percolating) cluster of sites joined by edges of weight , for there is a nonzero probability that the walker crosses the lattice without accumulating a single positive edge weight.
Regarding BGRW types A through D, Fig. 4 shows the critical points as function of bias strength for both, the bimodal disorder distribution Eq. (1) (see Fig. 4(a)) and the semicontinuous disorder Eq. (2) (see Fig. 4(b)). As evident from the figure, BGRWs of type A and C show a minimum at critical disorder parameter value at bias strengths in the range . In contrast to this, BGRWs of type B and D appear to reach a minimum only when the bias strength is zero. This observation holds for both, bimodal and semicontinuous disorder, although for the different disorder types the minima are located at different values and . For all the types of BGRWs considered, Tab. 1 lists the lowest value of , attained at the corresponding bias strength . This can be thought of as a measure of efficiently for the different walk dynamics. The lower the dynamical threshold , the more successful a walker upon traversing a graph aimed at finding a path of total negative weight (under the respective dynamics).
Another quantity that allows to characterize the different dynamics is the average number of steps taken by the different BGRW types until the walk terminates, see Fig. 5. This quantity can further be used to determine the fractal scaling dimension of the walks: given that a walk of steps spans an linear distance , its associated scaling dimension is defined via . Note that for the resulting paths, due to the possibility of loop erasure for some algorithms, the fractal dimension might be different. On a general basis, as soon as there is a bias that causes an effective drift of a walker (as in the BGRWs considered here), one can expect that on large enough lattices one trivially finds . Note that this is different when the environment itself is fractal, as for biased RWs on percolation clusters discussed in sect. II. Here, for completeness and to illustrate the precise dependence of the dynamics on the disorder parameter and different bias strengths , the rescaled number of steps taken during the lattice walks under the four BGRW dynamics are shown in Fig. 5 for the bimodal disorder distribution.
Considering type A BGRWs, see Figs. 5(a) and (b), it can be seen that the number of steps generally decreases when the density of negative edges is increased. This can be easily explained by taking into account the greedy behavior of the walker. In case that the walker only consumes the negative edges (replacing them with a standard positive value of 1; see rule (ii) discussed above), this behavior appears since, after each step taken, the probability of returning to the previous node is smaller than that of proceeding forward towards other directions that might have negative edges. In the case where a type A BGRW explores a graph with a low value of the disorder parameter and alters both, positive and negative edge weights (replacing them by edges with weight zero), an effective “repulsive effect” takes place because of the comparatively large probability that it might encounter excess edges with positive weights. Hence, the walker is likely to return to its previous location by traversing the edge with weight zero created while performing the last step. This repulsive effect is emphasized especially for a bias strength . As can be seen in Fig. 5(b), this repulsive effect leads to an increase in the number of steps taken by more than one order of magnitude in the low bias regime, compared to the case when only the negative bonds are altered. Also note that for different system sizes the data curves for the scaled number of steps fall on top of each other. Hence, the expected scaling dimension can be verified from the figures.
Regarding BGRWs of type B, illustrated in Fig. 5(c), the number of steps taken in the extreme cases (positive edge weights only) and (negative edge weights only), where there is no disorder at all, is the same, as one would expect. Regarding the bimodal disorder distribution, a minimum in the walk length is found when the disorder parameter assumes the value . This is caused by the fact that each node is, on average, endnode of one edge with negative edge weight that acts like a local trap for the walker. Only under the influence of the bias the walker can escape from such a trap. Until then, the walker has possibly traversed the edge with negative weight multiple times. After the loop erasure process (discussed for type B BGRWs in subsect. III.1), this results in an almost straightline walk. Increasing above also increases the likelihood of a node of having, on average, more than one incident edge with negative edge weight, thus preventing the walker from getting stuck (however, note that for the setup considered here, only above there is a system spanning cluster of negative edge weights). Although not shown here, in the case of a semicontinuous disorder distribution, the number of steps shows a similar behavior up to but instead of increasing afterwards it has a steady decrease until it reaches density which is caused by the lack of degeneracy in the bond weights.
A close similarity between BGRWs of type A and C can be seen upon comparison of the respective subfigures in Fig. 5. The most evident difference between the two dynamics is that type C BGRWs do not reach the minimum number of steps when the graph is fully equipped with negative edge weights. This effect appears due to the additional optimization overhead that takes advantage of the larger number of negative edge weights at large values of by increasing the length of the random walk.
For type D BGRWs, shown in Fig. 5(f), a behavior strikingly different from the other three BGRW dynamics can be observed. The respective lattice walks exhibit a minimum in the number of steps when no negative edges are present in the graph and attain a maximum when the graph is fully equipped with negative edge weights. This is a tell tale sign of an increase in the efficiency of traversing edges with a negative weight, while avoiding positive edge weights. Shorter lattice walks are performed at small values of and longer paths, possibly with many detours, are found at larger values of . The latter characteristic for the optimal paths in a graph predominantly equipped with negative edge weights.
It is interesting to note that the more efficient a random walk algorithm is at finding negative weighted paths, the higher the number of steps it will take at a low density of negative bonds if it alters both, negative and positive bonds, corresponding to rule (i) discussed above. As we observe, this effect is especially pronounced in the case of BGRWs of type D. The respective algorithm could not even be used efficiently for (no results for this case are shown).
As pointed out above, the expected scaling dimension for the BGRWs is . The number of steps taken during the lattice walks, normalized by the system size verifies this in the form for given values of and . Discrepancies in the curve overlap, corresponding to deviations from the expected scaling behavior, are accentuated at small values of , because the size of the walk gets closer to a dependent crossover length (similar to the crossover time discussed in sect. II). For walks much shorter than this crossover length the BGRWs exhibit a scaling dimension larger than . Only for walks longer than the crossover length, they assume the asymptotic scaling dimension . The simple explanation for this effect is that on short length scales the diffusive (albeit greedy) behavior dominates the configuration of the lattice walk, while on longer lengthscales the drift induced by the bias is dominant. Note that the scaling dimension of the lattice walks obtained here is very different from the fractal scaling dimension of the paths in the NWP problem. For the setup of the NWP problem, a fractal scaling dimension of was observed at Melchert and Hartmann (2008), increasing towards at (measured for the loopsonly setup) Melchert and Hartmann (2011).
iv.2 Ant colony optimization
In order to find values of and that reflect the underlying edgeweights, leading to an exploration in the vicinity of the “best” path found so far (referred to as intensification) but still allowing for an efficient exploration of many different paths (referred to as diversification), we performed several “calibration”runs of the ACO algorithm. During our simulation we used , , and . Below we report the results of a finitesize scaling analysis for the probability that the ACO algorithm finds a negatively weighted optimal path for different values of the fraction of negative edgeweights on square lattice graphs with . At each value of , we considered different realizations of the edge weights drawn from a bimodal disorder distribution (Eq. 1) in order to compute averages. The considered systems have periodic boundaries in the, say, vertical direction and open boundaries in the horizontal direction. Further, two extra “outer” nodes, i.e. the source node and the target node , are introduced. The source (target) node is connected to all nodes on the left (right) system boundary as shown in Fig. 1. Consequently, a path spans the system along the horizontal direction. Depending on the fraction of negative edge weights it might also wrap the system along the vertical direction. Fig. 6 illustrates exemplary paths found by the ACO algorithm for three different values of the disorder parameter on a square lattice of side length . In the figure, the black path indicates the best path found after the algorithm terminated and the linewidth of the other edges (colored grey) indicate how frequently the respective edge was visited by an agent. As evident from Fig. 6(a), corresponding to , if there are only few negative edge weights, the “best” path found by the ACO algorithm is rather straightlined. As the fraction of negative edge weights increases, see Figs. 6(b,c), corresponding to and , respectively, the paths roughen up and eventually also wind around the lattice in the vertical direction.
The characteristics of the weight of the best path found by the ACO algorithm as a function of the number of sweeps carried out by the algorithm is shown in Fig. 7(a). As can be seen in the curve corresponding to and after a short exploration time ( sweeps), the ACO algorithm maintains an average best pathweight . Since the underlying edgeweight distribution is bimodal, i.e. it allows for only, this reflects that at very small values of almost all edges contained in paths will have weight and the paths tend to have length . The situation changes upon increasing the value of . As can be seen from Fig. 7(a), already at and after a certain initial exploration time (), the ACO algorithm is able to identify paths with a negative weight so that .
Finally, Fig. 7(b) illustrates the scaling behavior of the probability , that the ACO algorithm returns a negative weighted path as a function of the fraction of negative edge weights. As can be seen from the inset of Fig. 7(b), the data curves for the different system sizes cross at . Below (above) and for increasing system size, the probability that the ACO algorithm returns a path weight tends to zero (one). As can be seen from the main plot of Fig. 7(b), can be rescaled using Eq. (5), similar to the paths found using the BGRW dynamics discussed previously. Considering the system sizes , a best data collapse (attained in the range on the rescaled axis) results in the estimates and with a quality of the data collapse Melchert (2009) (the numerical value of measures the mean–square distance of the data points to the master curve, described by the scaling function, in units of the standard error Houdayer and Hartmann (2004)). Performing the analysis for different intervals on the rescaled axis led to various estimates in the range . Further estimates of and can be obtained from the scaling behavior of the variance . It assumes a peak at an effective, dependent critical point at which (not shown). We find that seems to saturate at a value consistent with . In particular, we find (for the smaller system size we find ), obtained by fitting a Gaussshaped curve to the data of . The value of specifies the location of the peak of the fittingcurve and the errorbar is the respective fitting error. Further, the width of the Gaussshaped fitting function is consistent with , where we found and considering . Note that the value of the exponent found here is also in agreement with the one found for the BGRWs in the preceeding subsection.
V Conclusions
In the presented article we have investigated the principle characteristics of biased greedy random walks (BGRWs) on lattices with quenched disorder on the lattice edges. Four different types of BGRWs and an algorithm based on the ant colony optimization (ACO) heuristic were considered. Regarding the BGRWs, the precise configurations of the lattice walks constructed during the numerical simulations were influenced by two parameters: a disorder parameter that controls the amount of negative edge weights on the lattice and a bias strength that governs the drift of the walkers along a certain lattice direction. Focus of the presented study was the probability that, after termination, a lattice walk exhibits a negative weight as function of and . All four types of dynamics exhibit a phase transition with a characteristic disorder value above which lattice walks with a negative weight appear first. Further, the data curves for the probability that a walk with negative energy is found as function of could be rescaled well according to a simple scaling form. The respective scaling parameters in the extremal case of maximal bias , where the dynamics is effectively one dimensional, could be estimated from a simple renormalization argument. A lower bound on the value of can be obtained from the NWP problem in . Regarding those “static” simulations considering bimodal disorder on square lattices, negativeweighted paths appear first above the threshold with critical exponent . Considering the BGRW dynamics only, the type C BGRW (implementing rule (i) for the walker–edge weight interaction discussed in sect. III.1) yields the closest estimate, reading at bias strength , followed closely by type D lattice walks with at vanishing bias, see Tab. 1. Correspondingly, algorithm A (i) needs the largest number of steps to yield the best results. Interestingly, algorithm D is considerably faster, also compared to almost all variants.
The ACO based algorithm outperforms these estimates, yielding . Note that, in comparison to the type A and B BGRWs, which perform slightly worse, the dynamic rules governing the type C and D lattice walks and the ACO heuristic are rather involved. Also note that, albeit all the dynamic rules implement some degree of “local” optimization, they were not able to reach the threshold value obtained from the static simulations in terms of the NWP problem. Also the value of the critical exponent, for all types of BGRW and ACO dynamics, and the fractal behavior of the walks, hence of the paths, is different from the static negative weight paths. Thus static and dynamic behavior of the problem is different, as for many glassy systems.
Whether one can reach this static threshold and find the same critical exponents by means of simple dynamics that implement local optimization remains an intriguing open problem. This could be possible in principle, because the corresponding optimization problem can be solved in polynomial time. At least we could not succeed in doing so via the biased greedy random walks considered in the presented study.
Another interesting problem is whether a collective algorithm, which explores, e.g., the dynamics of an ensembles of paths connected to a heat bath at temperature , with the path weight being the energy, exhibits phase transitions as function of the disorder parameter as well. These transitions could be in the limit similar to the phase transitions found previously for the exact algorithms. For lower temperatures, maybe below a critical dynamical threshold, these transitions could be more similar to the transitions investigated in this work. Related studies are currently being performed.
Acknowledgements.
TLM acknowledges support from ESF through the project POSDRU 107/1.5/S/80765. OM acknowledges financial support from the DFG (Deutsche Forschungsgemeinschaft) under grant HA3169/31. The simulations were performed at the HPC Cluster HERO, located at the University of Oldenburg (Germany) and funded by the DFG through its Major Instrumentation Programme (INST 184/1081 FUGG) and the Ministry of Science and Culture (MWK) of the Lower Saxony State.References
 de Gennes (1971) P. G. de Gennes, La Recherche 7, 919 (1971).
 Straley (1980) J. P. Straley, J. Phys. C 13, 2991 (1980).
 Stauffer (1999) D. Stauffer, Physica A 266, 35 (1999).
 Melchert and Hartmann (2008) O. Melchert and A. K. Hartmann, New. J. Phys. 10, 043039 (2008).
 Apolo et al. (2009) L. Apolo, O. Melchert, and A. K. Hartmann, Phys. Rev. E 79, 031103 (2009).
 Melchert et al. (2010) O. Melchert, L. Apolo, and A. K. Hartmann, Phys. Rev. E 81, 051108 (2010).
 Melchert et al. (2011) O. Melchert, A. K. Hartmann, and M. Mézard, Phys. Rev. E 84, 041106 (2011).
 Norrenbrock et al. (2012) C. Norrenbrock, O. Melchert, and A. K. Hartmann (2012), preprint: arXiv:1205.1412.
 Claussen et al. (2012) G. Claussen, L. Apolo, O. Melchert, and A. K. Hartmann, Phys. Rev. E 86, 056708 (2012).
 Young (1998) A. P. Young, ed., Spin glasses and random fields (World Scientific, Singapore, 1998).
 Hartmann and Weigt (2005) A. K. Hartmann and M. Weigt, Phase Transitions in Combinatorial Optimization Problems (WileyVCH, Weinheim, 2005).
 Binder and Kob (2011) K. Binder and W. Kob, Glassy Materials and Diosordred Solids (World Scientific, Singapore, 2011).
 BenAvraham and Havlin (1982) D. BenAvraham and S. Havlin, J. Phys. A 15, L691 (1982).
 Gefen et al. (1983) Y. Gefen, A. Aharony, and S. Alexander, Phys. Rev. Lett. 50, 77 (1983).
 Pandey (1984a) R. B. Pandey, Phys. Rev. B 30, 489 (1984a).
 Argyrakis and Kopelman (1984) P. Argyrakis and R. Kopelman, J. Chem. Phys. 81, 1015 (1984), A summary of this article is available at papercore.org, see http://www.papercore.org/Argyrakis1984No2.
 White and Barma (1984) S. R. White and M. Barma, Journal of Physics A: Mathematical and General 17, 2995 (1984).
 Havlin and BenAvraham (1987) S. Havlin and D. BenAvraham, Advances in Physics 36, 695 (1987).
 Bouchaud and Georges (1990) J. P. Bouchaud and A. Georges, Phys. Rep. 195, 127 (1990).
 Argyrakis et al. (1984) P. Argyrakis, L. Anacker, and R. Kopelman, J. Stat. Phys. 36, 579 (1984), A summary of this article is available at papercore.org, see http://www.papercore.org/Argyrakis1984.
 Avramov et al. (1993) I. Avramov, A. Milchev, and P. Argyrakis, Phys. Rev. E 47, 2303 (1993), A summary of this article is available at papercore.org, see http://www.papercore.org/Avramov1993.
 Bunde et al. (1987) A. Bunde, H. Harder, S. Havlin, and H. E. Roman, J. Phys. A 20, L865 (1987).
 Arapaki et al. (1997) E. Arapaki, P. Argyrakis, I. Avramov, and A. Milchev, Phys. Rev. E 56, R29 (1997), A summary of this article is available at papercore.org, see http://www.papercore.org/Arapaki1997.
 Avramov et al. (1998) I. Avramov, A. Milchev, E. Arapaki, and P. Argyrakis, Phys. Rev. E 58, 2788 (1998).
 Dhar and Stauffer (1998) D. Dhar and D. Stauffer, Int. J. Mod. Phys. C 09, 349 (1998), A summary of this article is available at papercore.org, see http://www.papercore.org/Dhar1998.
 Kirsch (1999) A. Kirsch, Int. J. Mod. Phys. C 10, 753 (1999).
 Zhang et al. (2006) Z. X. Zhang, O. W.Z., X.W. Zou, and Z.Z. Jin, Physica A 360, 391 (2006).
 M. and Dhar (1983) B. M. and D. Dhar, J. Phys. C 16, 1451 (1983).
 Dhar (1984) D. Dhar, Journal of Physics A: Mathematical and General 17, L257 (1984).
 Pandey (1984b) R. B. Pandey, Phys. Rev. B 30, 489 (1984b).
 Skarpalezos et al. (2013) L. Skarpalezos, A. Kittas, P. Argyrakis, R. Cohen, and S. Havlin, Phys. Rev. E 88, 012817 (2013).
 Caracciolo et al. (1992) S. Caracciolo, A. Pelissetto, and A. D. Sokal, J. Stat. Phys. 67, 65 (1992).
 Kremer (1981) K. Kremer, Z. Phys. B 45, 149 (1981), A summary of this article is available at papercore.org, see http://www.papercore.org/Kremer1981.
 Grassberger (1993) P. Grassberger, J. Phys. A 26, 1023 (1993).
 Berretti and Sokal (1985) A. Berretti and A. Sokal, J. Stat. Phys. 40, 483 (1985).
 Smailer et al. (1993) I. Smailer, J. Machta, and S. Redner, Phys. Rev. E 47, 262 (1993).
 Majumdar (1992) S. N. Majumdar, Phys. Rev. Lett. 68, 2329 (1992).
 Grassberger (2009) P. Grassberger, J. Stat. Phys. 136, 399 (2009).
 Boyer et al. (2004) D. Boyer, O. Miramontes, G. RamosFernandez, J. L. Mateos, and G. Cocho, Physica A 342, 329 (2004).
 (40) For an overview of current developments regarding the ACO metaheuristic, see: http://www.acometaheuristic.org.
 Dorigo and Blum (2005) M. Dorigo and C. Blum, Theor. Comp. Sci. 344, 243 (2005).
 Dorigo (1997) M. Dorigo, IEEE Trans. Evol. Comp. 1, 53 (1997).
 Neumann and Witt (2010) F. Neumann and C. Witt, Theor. Comp. Sci. 411, 2406 (2010).
 Wilson (1971) K. G. Wilson, Phys. Rev. B 4, 3174 (1971).
 Wilson (1983) K. G. Wilson, Rev. Mod. Phys. 55, 583 (1983).
 Stanley (1999) H. E. Stanley, Rev. Mod. Phys. 71, S358 (1999).
 Melchert and Hartmann (2011) O. Melchert and A. K. Hartmann, Eur. Phys. J B 80, 155 (2011).
 Melchert (2009) O. Melchert, Preprint: arXiv:0910.5403v1 (2009).
 Houdayer and Hartmann (2004) J. Houdayer and A. K. Hartmann, Phys. Rev. B 70, 014418 (2004).