Ising Ferromagnets on Proximity Graphs with Varying Disorder of the Node Placement

We perform Monte Carlo simulations to determine the critical temperatures of Ising Ferromagnets (IFM) on different types of two-dimensional proximity graphs, in which the distribution of their underlying node sets has been changed systematically by means of a parameter σ. This allows us to interpolate between regular grids and proximity graphs based on complete random placement of nodes. Each edge of the planar proximity graphs carries a weighted ferromagnetic coupling. The coupling strengths are determined via the Euclidean distances between coupled spins. The simulations are carried out on graphs with N = 162 to N = 1282 nodes utilising the Wolff cluster algorithm and parallel tempering method in a wide temperature range around the critical point to measure the Binder cumulant in order to obtain the critical temperature for different values of σ. Interestingly, the critical temperatures depend partially non-monotonously on the disorder parameter σ, corresponding to a non-monotonous change of the graph structure. For completeness, we further verify using finite-size scaling methods that the IFM on proximity graphs is for all values of the disorder in the same universality class as the IFM on the two-dimensional square lattice.

The Ising model of ferromagnetism 1 is one of the most extensively studied models in statistical mechanics, because it features a continuous phase transition while being simple. The model consists of spins, which can be in one of two different states "+1" and "−1". The pair-wise interactions of spins, can be expressed by graphs and thus there is a constant interest in the critical behaviour of this model on different graph ensembles. The research provides the whole range from analytical results for regular lattices 2, 3 to numerical results for complex networks [4][5][6][7] . Further, there are many recent studies examining this model on different graphs [7][8][9][10][11] . This demonstrates the continuous interest in not only this model but also in related ones 12,13 . While most of the time the behaviour of the model is in the focus of the research, one can also use the model as a measuring probe to find properties of the underlying networks 14,15 .
An important type of graph is obtained from the Delaunay triangulation that finds application in, e.g., finite volume methods 16,17 . For a system of randomly placed Ising spins on a two-dimensional surface with neighbour relationships derived from a Delaunay triangulation, the critical temperature has been obtained and it has been confirmed that it is in the same universality class, i.e. it shows the same critical exponents as the two-dimensional Ising model on a square lattice [18][19][20] respectively in three dimensions on a cubic lattice 21 . This result is not trivial, since the Harris criterion suggests that random disorder is marginally important 22,23 for the two-dimensional Ising model.
The availability of results for this special case motivated us to extend and generalise the former research in two ways: First, we study two other types of irregular lattices, namely the Relative Neighborhood graph (RNG) 24 and Gabriel graph (GG) 25 , which are subgraphs of the Delaunay triangulation 26 and belong to the family of proximity graphs. The objective of proximity graphs is to connect nodes which are spatially close to each other, hence they are suited to generalise problems defined on regular lattices with nearest neighbour relationships. This family of graphs finds application in geographic variation studies in biology [27][28][29] , as potential candidates for "virtual backbones" in ad-hoc networks [30][31][32][33][34] and in machine learning and pattern recognition 35 . Also, their statistical properties have been under scrutiny recently 36,37 .
As a second generalisation and extension, we consider an ensemble of node sets depending on a tunable parameter σ. By changing this parameter we can interpolate the configurations between a square lattice and nodes distributed by a Poisson point process. A similar approach has recently been used 38 to introduce an ensemble of Travelling Salesperson Problems (TSP) which interpolates between an easy to solve circular arrangement and a hard to solve random placement of the nodes. The increase of available computing power in recent years made it possible to perform the simulations at similar number of sites but increased number of realisations of the randomness and many values for σ in comparison to previous studies for the large σ limit 19,20 .
For any node set the construction rule of a proximity graph defines the corresponding edge set. Thus, the parameter σ influences the regularity of the resulting proximity graph. Below, we measure the critical temperatures for a range of σ. While mean field theory T c = JK and the Bethe approximation , which is practically linear for > K 3  , both suggest a linear dependence between T c and the mean number of neighbours K, we find approximately a power-law relationship with an exponent which is small but certainly larger than one. Additionally, also to confirm the quality of our data, we verify by applying a common finite-size scaling analysis to obtain the critical exponents and by considering the two-point finite-size correlation function that the universality class does not depend on the graph type and does not change when varying σ.
The remainder of this work is organised as follows: First the RNG, GG and the Ising model are introduced. Then the details and results of our simulations are discussed. At the end a conclusion of this article is given.

Model Proximity
Graphs. An undirected graph G(V,E) consists of a set of nodes V and edges E⊂V (2) . Two nodes connected by an edge are called neighbours. Proximity graphs are defined in a metric space. In this article a two-dimensional space with Euclidean metric is chosen, because it is the most common case and easy to visualise. Though every metric in any dimension can be used, in principle.
One of the proximity graphs we consider here is the Relative Neighbourhood graph (RNG) 24 . Within this graph two nodes i and j with distance d ij will be connected, if no other node is located in a well defined area called lune. The lune is defined as the intersection area of two disks with radius r = d ij , whose centres are on node i and j, respectively. This means that the edge {i, j} will be part of the RNG, if the condition ij ik kj is fulfilled. For the sake of clarity, the lune of two nodes in regard to the RNG is sketched in Fig. 1a (hatched region). An example of the RNG is given in Fig. 1b. For a Poisson point process in a square, i.e. nodes placed on uniformly distributed, independent coordinates, its mean degree is K = 2.5576(3) 36 , i.e. the mean number of neighbours of each node. The second proximity graph we study is the Gabriel graph (GG) 25 . In the GG two nodes i and j with distance d ij will be connected, if the disk which has the connecting line between i and j as its diameter contains no further nodes k∈{i,j}, i.e. the condition is fulfilled. The lune is illustrated in Fig. 1a (cross hatched region). An example of the GG is given in Fig. 1c. Note that the lune of the GG is completely enclosed in the lune of the RNG, therefore the RNG is a subgraph of the GG and, as a consequence, all edges belonging to the RNG also belong to the GG. For a Poisson point process in a square, its mean degree is K = 4 37 -the same as the square lattice. Note that by construction RNG and GG have no crossing edges, i.e. they are planar. Furthermore they are connected, i.e. consist only of one connected component. For any given set of nodes located in a two-dimensional plane these graphs can be constructed by an algorithm 36 featuring a time complexity of  N N ( log ). In this algorithm the Delaunay triangulation, which can be  40 and which is a supergraph of the RNG and GG, is obtained first. Subsequently its edges, of which there are = E N ( )  , are checked in regard to the aforementioned construction rules. A simple construction of proximity graphs can be based on a uniform and random placement of nodes. To allow for a more complex behaviour, we apply a more general approach. Here, first the N = L 2 nodes are placed to form a regular square lattice with periodic boundary conditions. Next, the nodes are displaced, introducing geometric disorder resulting in a non-regular graph structure. The horizontal displacement Δx and vertical displacement Δy of each node is drawn from a Gaussian distribution with zero mean and standard deviation σ. Therefore, via the parameter σ we are able to interpolate between fully regular and fully random placement. The influence of σ on the structure of a proximity graph has two main effects which are clearly visible in Fig. 2 for the RNG. First, the number of edges changes as will be examined in more detail later on. Second, the variance of the length of the edges increases, i.e. there are some quite long and some very short edges in configurations with increasing σ. The Ising Model. Commonly, the IFM is studied on a square lattice with lateral length L and N = L 2 nodes with periodic boundary conditions. Each node i has an Ising spin s i ∈{−1, +1} and interacts with its neighbours according to the Hamiltonian 〈i, j〉 refers to nodes i and j which are nearest neighbours, i.e. they are connected by an edge. The coupling strength between i and j is given by J ij .
In the generalised graph we use here, the interaction structure is not regular. To account for the varying distances, we choose the coupling constants J ij , which depend on the geometric distance between the connected pair of nodes. More precisely, the coupling constant of an edge {i, j} is where d ij is the Euclidean distance between node i and j. Following the work of Lima et al. 20 . the free parameter α is set to α = 0.5. Note that distances smaller and larger than one can appear. For σ = 0 all d ij = 1 and therefore all J ij = 1 for every edge {i, j} ∈ E.

Results
The

Reproduction of the Critical Exponents. First, we want to confirm that the Ising model on both planar
proximity graphs belong to the same universality class as the Ising model on a square lattice. This can be since dν = 2 for the Ising model and therefore, according to the Harris criterion 22 , the disorder is not relevant, or at least border-line, i.e. its influence is hard to detect 46 . Therefore not only the critical exponents were obtained, but also a sensitive test using the two-point finite-size correlation function was done. While this examination will not yield any surprising results, it demonstrates the reliability of the simulation data and ensures the correctness of the simulations. The critical exponents were determined by a finite-size scaling (FSS) analysis 47 . Therefore we looked at the magnetisation m, susceptibility χ and Binder cumulant g, which are defined as Here, 〈…〉 denotes the thermal average and […] avg the average over the disorder realisations. For all quantities, we obtained statistical errors via bootstrap resampling 48,49 . We determine the exponents by using the finite-size scaling assumption 50, 51 for phase transitions of 2nd order, e.g., for the susceptibility χ: where ⋅  C( ) is some unknown scaling function. That means, for the right choice of ν, γ and T c and appropriately scaled axes, the measurements at different system sizes L collapse on the scaling function as shown in Fig. 3a. A similar formula is valid for m (containing β) and g (which is dimensionless, i.e. has no prefactor in front of the scaling function). Values for the exponents and T c are obtained in an automated and reproducible way 52 . They are visualised in Table 1 and are within errorbars consistent with the exactly known values for the square lattice IFM.
A more sensitive method 53,54 to test whether two models are in the same universality class includes the analysis of the two-point finite-size correlation function Figure 3. (a) Example of a data collapse of the susceptibility χ in regards to the RNG at σ = 1 to get the corresponding critical exponent γ. The same procedure is used on the Binder cumulant g and the magnetisation |m| to determine the critical exponents ν and β and the critical temperature T c (see Table 1). The inset shows the same data with unscaled axes. (b) Binder cumulant g as a function of the two-point finite-size correlation function divided by system size L. The data points for different system sizes are drawn with different symbols, while the data points for different graph ensembles are drawn with different colours. All points fall on one curve within errorbars.  which is in this case only evaluated over one direction, thus x j is the position of the node in the horizontal direction.
We have plotted the Binder cumulant g as a function of ξ/L for different values of σ, L and T as shown in Fig. 3b. We can observe that within error bars all data points fall on the same curve, which confirms that the universality class does not change by varying σ.
Behaviour of the Critical Temperature. This study examines the critical temperatures T c for different values of σ. In order to locate T c we considered the Binder cumulants for different system sizes L, because they are intersecting at the critical point T c 55 as shown in Fig. 4a. Since g is only measured for discrete values of T, the points have to be interpolated to find the intersection. Therefore a cubic spline 56 interpolation is calculated for the measured points. An even better method would be the multiple histogram method 45,57 . But the simple cubic spline interpolation seems to deliver results that are good enough, since the comparison for the analytically known case σ = 0, i.e. square lattice, with T c = 2.2691… 2 suggests that the cubic splines offer good accuracy, since both = . T (0) 2 2690 (2) c RNG and = . T (0) 2 2690 (3) c GG are perfectly compatible with the reference value and offer reasonable precision.
In Fig. 4b the critical temperatures are shown as a function of the disorder strength. Interestingly, we first observe that for the GG there occurs a jump from T c (0) = 2.2690(3) to a value close to the measured value T c (0.03) = 2.8511 (7). This is easily explainable, because the square lattice at σ = 0 is a very special case for the GG. At very small deviations from the square lattice new edges arise in the GG, this is visualised in Fig. 5a-c. These new edges lead to a stronger coupling of spins such that the system stays ferromagnetic at higher temperatures. After an initial increase of the GG's critical temperature with growing disorder, a maximum is reached and the critical temperature drops. Thus, the dependence is non-monotonous. This can be understood when analysing the dependence of the average number K of neighbours as a function of the disorder strength σ, as shown in Fig. 4b. The behaviour of T c and K is very similar: K(σ) exhibits also a jump at σ = 0 as well as a peak close to σ = 0.2. This is confirmed when plotting T c as a function of K, see Fig. 4c. While the Bethe approximation predicts an almost linear relationship T c ∝ K, both proximity graphs show a power-law relationship T c = aK b between the degree K and the critical temperature T c with an exponent b ≈ 1.3. The same analysis is applied to a further simulation with all edge weights set to J = 1, i.e. α = 0, which also shows b ≈ 1.3. This can be interpreted as the edges chosen by the proximity graphs are more efficient to couple the spins than in the Bethe lattice. This stronger coupling is expressed by a higher T c . Note that the σ = 0 point is not part of this power law, which is expected, since it is a special corner case as mentioned above. Mind however, that the power-law relation is only approximately correct, since especially the RNG shows deviations from this form for low degrees K.
For large values of σ the node configurations approach the limit of a Poisson point process, hence T c approaches a limit value, which was obtained from a fit as T c (σ = ∞) = 2.10(3) for the GG and T c (σ = ∞) = 1. 19 RNG. A further interesting feature which highlights the importance of the structure of the graph for the critical temperature is the fact that the square lattice has a slightly different critical temperature than the GG in the Poisson point process limit although they share the same average degree K = 4 37 . For the RNG, the critical temperature T c of the RNG decreases smoothly and monotonously with increasing value σ. Therefore, the RNG exhibits neither a jump at σ = 0 nor a peak. This is reflected by the behaviour of K(σ) originating from a decrease of the number of edges, see Fig. 2 and Fig. 4b. Note also that for the RNG the power law between critical temperature and number of neighbours is well visible and also shows also an exponent of b ≈ 1.3 except for the smallest K. This is plotted in the inset of Fig. 4c.

Conclusion
We have performed Monte Carlo simulations of an Ising Ferromagnet in a two-dimensional Euclidean plane with neighbour relationships defined by two different proximity graphs, the Relative Neighbourhood graph and Gabriel graph. Controlled by a disorder parameter σ we can smoothly interpolate between a regular square lattice node arrangement and a completely random node arrangement.
By means of finite-size scaling analysis and by studying the two-point finite-size correlation function, we confirmed the expectation that the universality class does not change for different graph types and for different values of σ. Furthermore, we have monitored the critical temperature T c as a function of σ. For GG we observe a jump and a non-monotonous behaviour, reflecting non-trivial changes of the graph topology. On the other hand, for RNG, T c stays approximately constant for small values of σ and decreases to some limit value for  large disorder. In fact, T c strongly depends on the degree K. It turns out that there is a power-law relationship between T c and K for both proximity graphs. Even if the degree K = 4 of the GG at σ = 0 is equal to the degree for large values of σ, the critical temperature is lower at this point. Therefore, it becomes obvious that the average number of interaction partners is the dominating but not the sole influence for the location of the critical point.