High-Temperature Series Expansions for Random Potts Models

We discuss recently generated high-temperature series expansions for the free energy and the susceptibility of random-bond q-state Potts models on hypercubic lattices. Using the star-graph expansion technique quenched disorder averages can be calculated exactly for arbitrary uncorrelated coupling distributions while keeping the disorder strength p as well as the dimension d as symbolic parameters. We present analyses of the new series for the susceptibility of the Ising (q=2) and 4-state Potts model in three dimensions up to order 19 and 18, respectively, and compare our findings with results from field-theoretical renormalization group studies and Monte Carlo simulations.


Introduction
Systematic series expansions [1] for statistical physics models defined on a lattice provide an useful complement to field-theoretical renormalization group studies and large-scale numerical Monte Carlo simulations. This is in particular true when studying phase transitions and critical phenomena of quenched, disordered systems. In the field-theoretic treatment [2] the necessary average over disorder realizations at the level of the free energy requires the application of the so-called "replica trick" which loosely speaking introduces n different, interacting copies of the original system, with the formal limit n → 0 taken at the end. In the numerical approach the average over a large but finite number of different disorder realizations can, at least in principle, be performed explicitly but is very time consuming such that only few points in the vast parameter space of the systems can be sampled with realistic effort. Moreover extrapolations of the data on finite lattices to the infinite-volume limit are required. Using high-temperature series expansions, on the other hand, one can obtain for many quantities exact results up to a certain order in the inverse temperature. Here the quenched disorder is treated exactly and the infinite-volume limit is implicitly implied. Moreover, one can keep the disorder strength p as well as the dimension d as symbolic parameters and therefore analyse large regions of the parameter space of disordered systems. The critical part of the series expansion approach lies in the extrapolation techniques which are used in order to obtain information on the phase transition behaviour from the finite number of known coefficients of the high-temperature series. While for pure systems this usually works quite well, one can question the use of these extrapolation techniques in disordered systems, where the singularity structure of the free energy or susceptibility may be very complicated, involving Griffiths-type singularities or logarithmic corrections [3].
Pure Potts models show either first-or second-order phase transitions, depending on the dimension d and the number of states q. Since in the second-order case the specific-heat exponent α is non-negative for this class of models, the Harris criterion [4] suggests for the corresponding disordered systems either the appearance of a new random fixed point (d = 2, q = 3, 4 and d = 3, q = 2) or logarithmic corrections to the pure fixed point (d = 2, q = 2). At first-order phase transitions, randomness softens the transitions [5]. For d = 2 even infinitesimal disorder induces a continuous transition [6,7], whereas for d = 3, q > 2 a tricritical point at a finite disorder strength is expected [8].
In this work we studied these scenarios by means of "star-graph" high-temperature series expansions where the disorder average can be taken at the level of individual graphs. Using optimized cluster algorithms for the symbolic, exact calculation of spin-spin correlators on finite graphs with arbitrary inhomogeneous couplings, we obtained series expansions for the free energy and susceptibility in the inverse temperature up to order 19 respectively 18 for bond-diluted Ising and Potts models in dimensions d ≤ 5, and up to order 17 in arbitrary dimensions. Here we shall focus on analyses of these series in three dimensions where a direct comparison with field-theoretic renormalization group studies and recent Monte Carlo simulations is possible.

Model
The ferromagnetic disordered q-state Potts model on hypercubic lattices Z d is defined by the partition function where β = 1/k B T is the inverse temperature, J ij are quenched (non-negative) nearest-neighbour coupling constants, the spins can take the values s i = 1, . . . , q, and δ .,. is the Kronecker symbol. In our series expansion the combination will be the relevant expansion parameter which in the Ising case (q = 2) simplifies to v ij = tanh(βJ ij /2). In the symmetric high-temperature phase, the susceptibility associated with the coupling i h i (qδ s i ,1 − 1)/(q − 1) to an external field h i is given for a graph with N spins by summing over all two-point correlations, Here the brackets [. . .] av indicate the quenched disorder average which in our case is taken over an uncorrelated bimodal distribution of the form Besides bond dilution (R = 0), which will be in the focus of the present work, this also includes random-bond ferromagnets (0 < R < 1) and the physically very different class of spin glasses (R = −1) as special cases. Other distributions such as Gaussian distributions can, in principle, also be considered with our method.

Series generation methodology
In this section we briefly review the main technical ingredients necessary for our high-temperature series study. We begin with a few basic notations from graph theory. A graph of order E consists of E links connecting N vertices. We consider only connected, undirected graphs that are simple: no link starts and ends at the same vertex (no tadpoles) and two vertices are never connected by more than one link. Subgraphs are defined by the deletion of links. In this process, isolated vertices can be dropped. Since each link may be present or absent, a graph of order E has 2 E (not necessarily non-isomorphic) subgraphs. These subgraphs may consist of several connected components and are called clusters. If the deletion of one vertex renders the graph disconnected, such a vertex is termed articulation point. The "star graphs" we are considering here are thus just defined by the absence of such articulation points. A graph is bipartite if the vertices can be separated into red and black vertices so that no link connects two vertices of the same color. Equivalently, all closed paths in the graph consist of an even number of links. Hypercubic lattices are evident examples for bipartite graphs.
There are a couple of well-established methods [1] known for the systematic generation of high-temperature series expansions which differ in the way relevant subgraphs are selected or grouped together. A recently developed alternative method [9] exploits ideas from so-called finite-lattice methods, usually employed before for the generation of low -temperature series. Using a clever reformulation of the method, Arisue et al. [10] succeeded to generate a very impressive 32th order world-record high-temperature susceptibility series for the pure Ising model in three dimensions.
For the class of classical O(N) spin models without disorder, quite long series (up to order β 25 ) have also been produced by linked-cluster expansions [11]. This technique also allows one to obtain series for more involved observables (such as the second moment of the spin-spin correlation function yielding the correlation length) which have no star-graph expansion. Furthermore, it works with free embeddings of graphs into the lattice which can be counted orders of magnitude faster than the weak embedding numbers needed by the star-graph technique. Nonetheless, the linked-cluster method has not yet been applied to problems with quenched disorder.
The star-graph method can be adopted to systems involving quenched disorder [12,13] (as also can the no-free-end method [14]) since it allows one to take the disorder average on the level of individual graphs. The basic idea is to assemble the value of some extensive thermodynamic quantity F on a large or even infinite graph from its values on subgraphs: Graphs constitute a partially ordered set under the "subgraph" relation. Therefore, for every function F (G) defined on the set of graphs exists another function W F (G) such that F (G) = g⊆G W F (g), for all graphs G. This function can be calculated recursively via where (G : Z d ) denotes the weak embedding number of the graph G in the given lattice structure [15].
The following observation makes this a useful method: Let G be a graph with an articulation vertex where two star subgraphs G 1,2 are glued together. Then ). An observable F for which this property is true on arbitrary graphs with articulation points allows a star-graph expansion. All non-star graphs have zero weight W F in the sum for F (Z d ). It is easy to see that the (properly normalized) free energy log Z has this property and it can be proved [13] that the inverse susceptibility 1/χ has it, too, even for arbitrary inhomogeneous couplings J ij . This restricts the summation for F (Z d ) to a sum over star graphs. The linearity of the recursion relations then enables the calculation of quenched averages over the coupling distribution on the level of individual graphs.
The resulting recipe for the susceptibility series is: • Graph generation and embedding number counting.
• Calculation of Z(G) and the correlation matrix for all graphs as polynomials in E variables v ij defined in (2).
• Inversion of the Z polynomial as a series up to the desired order.
• Averaging over quenched disorder, resulting in a matrix of polynomials in (p, v). • Inversion of the matrix N nm and subgraph subtraction, • Collecting the results from all graphs, Algorithmically the most cumbersome part of this recipe is the first step, the generation of star graphs and calculation of their (weak) embedding numbers. The graph generation is usually done by recursively adding nodes and edges to a list of smaller graphs. To make sure that no double counting occurs this requires an isomorphism test, i.e., the decision whether two given adjacency lists or adjacency matrices describe the same graph modulo relabelling and reordering of edges and nodes. We employed the NAUTY package by McKay [16] which allows very fast isomorphism tests by calculating a canonical representation of the automorphism group of the graphs. By this means, we classified for the first time all star graphs up to order 19 that can be embedded in hypercubic lattices, see Table 1. As with any series expansion, the effort grows exponentially with the maximal order of the expansion, rendering each new order roughly as "expensive" as all previous orders taken together. This is illustrated in Fig. 1 where already the number of star graphs is seen to grow exponentially as a function of the links E. The exponential fit in the range E = 13 -19 suggests that the number of star graphs increases roughly by a factor of 2.8 in each of the next higher orders, predicting about 65 000 different star graphs with E = 20 and about 180 000 with E = 21.
For each of these graphs we calculated their (weak) embedding numbers for ddimensional hypercubic lattices (up to order 17 for arbitrary d, order 18 (general q-state Potts) and 19 (Ising) for dimensions d ≤ 5). Two typical results are depicted in Fig. 2. For the embedding count we implemented a refined version of the backtracing algorithm by Martin [15], making use of a couple of simplifications for bipartite hypercubic lattices Z d . After extensive tests to find the optimal algorithm for the "innermost" loop, the test for collisions in the embedding, we ended up using optimized hash tables.
The second step of the series generation requires the exact calculation of the partition function and the matrix of correlations M nm for each star graph with arbitrary symbolic couplings J ij defined on the E ≤ 19 edges. The crucial observation is that this can be done most efficiently by using the cluster representation where the sum goes over all clusters C ⊆ G, e is the number of links of the cluster and c the number of connected components of C. The reduced partition function Z ≡ Zq E−N / ij (e βJ ij −1+q) is normalized such that log Z has a star-graph expansion. Similarly, the calculation of the susceptibility involves the matrix of correlations where the sum is restricted to all clusters C nm ⊆ G in which the vertices n and m are connected. This representation essentially reduces the summation over q N states to a sum over 2 E clusters which, compared with previous implementations, results is a huge saving factor in computing time (of the order of 10 6 ). Further improvements result if the 2 E clusters belonging to a graph are enumerated by Gray codes [17] such that two consecutive clusters in the sum (5) differ by exactly one (added or deleted) link. In the Ising case q = 2 another huge simplification takes place since only clusters where all vertices are of even degree contribute to the cluster sum. Since general purpose software for symbolic manipulations turned out to be too slow for our purposes, we developed a C ++ template library using an expanded degree-sparse representation of polynomials and series in many variables. For arbitrary-precision arithmetics the open source library GMP was used. Finally, for the case of bond dilution (R = 0 in (4)) considered here, we made use of the fact that the disorder average is most easily calculated via 4. Series analysis: techniques and results

Bond-diluted 3D Ising model
Disordered magnetic systems belonging to the 3D Ising model universality class have been studied extensively in experiments [18][19][20] and also by field theoretical and numerical methods. A comprehensive compilation of recent results can be found in [21], showing a wide scatter in the critical exponents of different groups, presumably due to large crossover effects.
Our high-temperature series expansion for the susceptibility up to order 19 is given with coefficients as polynomials in p, χ(v) = n a n (p)v n [22]. Therefore it should be well-suited for the method of partial differential approximants [23] which was successfully used to analyse series with an anisotropy parameter describing the crossover between 3D Ising, XY and Heisenberg behaviour [24]. But this method was unable to give conclusive results. Therefore we confined ourselves to a singleparameter series for selected values of p.
The ratio method assumes that the expected singularity of the form  is the closest to the origin. Then the consecutive ratios of series coefficients behave asymptotically as r n = a n a n−1  [25] (where T c goes to 0, v c to 1) the series is clearly ill-behaved, related to the exp(1/T ) singularity expected there. Besides that, the slope (related to γ) is increasing with p. The widely used DLog-Padé method consists in calculating Padé approximants to the logarithmic derivative of χ(v), The smallest real pole of the approximant is an estimation of v c and its residue gives γ. The results presented in Table 2 are the averages of 45 -55 different Padé approximants for each value of p, with the error in parentheses indicating the standard deviation. The scattering of the Padé approximants increases with p, getting again inconclusive near the percolation threshold. Nevertheless, up to about p = 0.6 the series estimates for v c respectively T c are in perfect agreement 1 with the Monte Carlo (MC) results of Ref. [26]. This is demonstrated in Fig. 4 where also the (properly normalized) mean-field and effective-medium approximation [27] are shown for comparison.
The critical exponent γ, as provided by this method, apparently varies with the disorder strength. More sophisticated analysis methods, such as inhomogeneous  (60) differential approximants [28,29], the Baker-Hunter method [30] or the methods M1 and M2 [31], especially tailored to deal with confluent singularities as one would expect in a crossover situation, give improved results in the pure (p = 0) case but do not essentially change the results in the presence of disorder. Thus, while for theoretical reasons we still find it likely that the variation of γ with the disorder strength can be attributed to neglected or insufficiently treated correction terms, it proved clearly impossible to verify this effect in the series analysis. In fact, a plot of γ vs. p does not even show an indication of a plateau. In the   [26]. For comparison also the (properly normalized) mean-field and effective-medium approximations are shown.

Bond-diluted 4-state Potts model
In three dimensions the 4-state Potts model exhibits in the pure case a strong first-order transition [36] which is expected to stay first order up to some finite disorder strength, before it gets softened to a second-order transition governed by a disorder fixed point.
In the latter regime we are interested in locating power-law divergences of the form (8) from our susceptibility series up to order 18 [37,38]. To localize a first-order transition point, however, a high-temperature series alone is not sufficient since there the correlation length remains finite and no critical singularity occurs. In analysing series by ratio, Padé or differential approximants, the approximant will provide an analytic continuation of the thermodynamic quantities beyond the transition point into a metastable region on a pseudo-spinodal line with a singularity T * c < T c and effective "critical exponents" at T * c . Again we first employed the ratio method which is the least sophisticated method of series analysis, but usually it is quite robust and gives a good first estimate of the series behaviour. Figure 5 shows these ratios for different values of p. They behave qualitatively similar to the Ising model case (oscillations caused by the antiferromagnetic singularity at −v c , strong influence of the percolation point at p c ≈ 0.75). Notice that the slope (∝ γ − 1) is increasing with p, changing from γ < 1 to γ > 1 around p = 0.5. Figure 6 compares the critical temperature, estimated from an average of 25 -  30 Padé approximants for each value of p, 2 with the results of recent Monte Carlo simulations [39]. For small p, in the first-order region, the series underestimates the critical temperature. As explained above, this is an estimate not of T c but of T * c . Between p = 0.3 and p = 0.5, the estimates confirm, within errors, the Monte Carlo results, indicating that now both methods see the same second-order transition. Beyond p = 0.5, the scatter of different Padé approximants increases rapidly, related to the crossover to the percolation point.
The situation is more complicated with respect to the critical exponent γ. The DLog-Padé analysis gives inconclusive results due to a large scattering between different Padé approximants, as shown in Fig. 7. One possible reason for this failure is the existence of confluent singularities. The dots in Eq. (8) indicate correction terms which can be parametrized as follows: where ∆ i are the confluent correction exponents. Among the various sophisticated analysis methods (inhomogeneous differential approximants [28,29] and the methods M1 and M2 [31]), in the case at hand, the Baker-Hunter method [30] appeared to be the most successful, giving consistent results at larger dilutions p > 0. 35 where the leading-term DLog-Padé analysis failed. The Baker-Hunter method assumes that the function under investigation has confluent singularities which can be transformed into an auxiliary function g(t) that is meromorphic and therefore suitable for Padé approximation. After the substitution z = z c (1 − e −t ) we expand F (z(t)) = n c n t n and construct the new series such that Padé approximants to g(t) exhibit poles at t = 1/λ i with residues −A i /λ i . This method is applied by plotting these poles and residues for different Padé approximants to g(t) as functions of z c . The optimal set of values for the parameters is determined visually from the best clustering of different Padé approximants, as demonstrated in Fig. 8. Using this method, our results for the critical exponent γ are plotted in Fig. 9. They show an effective exponent monotonically increasing with p but reaching a plateau at γ = 1 for dilutions between p = 0.42 and p = 0.46. The following sharp increase is to be interpreted as due to the crossover to the percolation fixed point at p c ≈ 0.75, T c = 0, where a χ ∼ exp(1/T ) behaviour is expected. It is well known (see, e.g., Ref. [40]) that series analysis in crossover situations is extremely difficult. If the parameter p interpolates between regions governed by different fixed points, the exponent obtained from a finite number of terms of a series expansion must cross somehow between its universal values, and does this usually quite slowly. Therefore it does not come as a surprise that the Monte Carlo simulations quoted above see the onset of a second order phase transition already for smaller values of the disorder strength p. The mere existence of a plateau in γ eff (p), however, is an indication that here truly critical behaviour is seen. It is governed by a fixed point for which we obtain γ = 1.00 (3). Here, as always in series analyses, the error estimates the scattering of different approximants.

Discussion
We have implemented a comprehensive toolbox for generating and enumerating star graphs as required for high-temperature series expansions of quenched, disordered systems. Monte Carlo simulations of systems with quenched disorder require an enormous amount of computing time because many realizations have to be simulated for the quenched average. For this reason it is hardly possible to scan a whole parameter range. Using high-temperature series expansions, on the other hand, one can obtain this average exactly. Since the relevant parameters (degree of disorder p, spatial dimension d, number of states q, etc.) can be kept as symbolic variables, the number of potential applications is very large.
Here we presented analyses of the susceptibility series for the three-dimensional bond-diluted Ising and 4-state Potts model. The resulting phase diagrams in the p-T -plane are in very good agreement with recent Monte Carlo results. As far as the critical exponent γ is concerned, however, large crossover effects render a reliable determination from series expansions up to order 19 respectively 18 very difficult.
In the Ising case we estimate values that are clearly different from the pure case but exhibit a pronounced dependence on the degree of dilution. For the 4-state Potts model with its strong first-order phase transition in the pure case, the singularity structure of the disordered model is even more involved. Still, by comparing the series expansions with numerical data we can identify signals for the onset of a softening to a second-order transition at a finite disorder strength.