Two-Dimensional Quantum Gravity – A Laboratory for Fluctuating Graphs and Quenched Connectivity Disorder

This paper gives a brief introduction to using two-dimensional discrete and Euclidean quantum gravity approaches as a laboratory for studying the properties of fluctuating and frozen random graphs in interaction with “matter fields” represented by simple spin or vertex models. Due to the existence of numerous exact analytical results and predictions for comparison with simulational work, this is an interesting and useful enterprise.


Introduction
Field theoretical formulations of Einstein gravity are known to be perturbatively non-renormalizable.Over the past decades, this has prompted active research into several constructive approaches to non-perturbative quantization prescriptions [1].Among them, Regge calculus [2] and the dynamical triangulations model in its Euclidean [3] and more recently Lorentzian [4] versions have been extensively studied over the last 20 years.The basic idea of both approaches is the same as in Feynman's formulation of quantum mechanics in terms of path integrals [5].While in quantum mechanics one sums over all paths a particle can take from a point x 0 at time t 0 to a point x 1 at time t 1 to compute the corresponding probability amplitude, in the gravity context one describes the quantum fluctuations of space-time by performing a functional integral over an ensemble of discrete, simplicial manifolds [6].In Regge calculus [2] the connectivities of the discretised piecewise linear manifolds are fixed and the edge lengths are the dynamical degrees of freedom.In the dynamical triangulations model, on the other hand, the situation is reversed: here the edge lengths are kept fixed, but now the connectivities are allowed to vary dynamically from vertex to vertex [6].This latter case allows for exact solutions.
From the viewpoint of any additional spin models coupled to the triangulations, the fluctuating manifolds act as a special kind of annealed disorder.By freezing in randomly selected manifolds and thus turning off the back-reaction of the matter fields, this class of systems leads quite naturally to the problem of spin models subject to quenched geometric disorder.In the Regge case with varying edge lengths, this corresponds to a specific type of quenched random-bond disorder.For the dynamical triangulations model, on the other hand, quenched connectivity disorder can be studied in this way.Based on the exact solutions for the annealed case, theoretical conjectures have also been made for the quenched situation.The complementarity of analytical and numerical methods is one of the main merits of the dynamical triangulations approach.
After a brief introduction to the two formulations of two-dimensional Euclidean quantum gravity, this paper will focus on the statistical physics interpretation of spin and vertex models coupled to fluctuating or quenched quantum gravity graphs.Both analytical and numerical results will be discussed and compared with each other.

Regge calculus
The Regge formulation of quantum gravity [2,7] stays relatively close to the continuum formulation, which for instance for two-dimensional (2D) Euclidean R 2 -gravity would be defined through the partition function with the gravitational action taken as where λ is the cosmological constant, g is the metric tensor and R is the scalar curvature.The functional integral measure Dµ(g) controlling the fluctuations of the manifolds described by the metric tensor g is usually taken as the DeWitt measure [8].This formally implements the basic idea for going from the classical to the quantum world à la Feynman: instead of finding a classical solution by optimising the action by the variational principle, one integrates over all possible manifolds parametrised by the metric tensor g, in close analogy to the sum-over-paths prescription in ordinary particle quantum mechanics.
To render Z(A) in (1) computable in practise, some discretisation is necessary.Regge's discretisation program [2] consists of replacing a given continuum manifold by piecewise linear manifolds, whose internal geometry is flat.This procedure works for any space-time dimension and for metrics of arbitrary signature.Originally it was applied as a computational tool to the classical optimisation problem.Here we restrict ourselves in the quantum context to the simplest case of two dimensions and Euclidean signature.Typically (but not necessarily) one considers piecewise linear manifolds with fixed connectivities.The dynamic degrees of freedom are then the edge lengths of the simplicial discretisation.
In two dimensions this procedure is most easily visualised by choosing a triangulation of the surface under consideration, where each triangle then represents a building block of the piecewise linear manifold.The net of triangles is itself a two-geometry, with singular (non-differentiable) points located at the vertices of the net, where several triangles meet.A vector that is linearly transported around these vertices experiences in the presence of curvature a rotation by the deficit angle δ i = 2π − t⊃i θ i (t), where θ i (t) is the dihedral angle at vertex i.For the area assignment one usually uses a barycentric decomposition, where A i = t⊃i A t /3 denotes the barycentric area with A t being the area of a triangle t.In order to derive the transcription from the continuum to the discretised formulation, one identifies the following continuum quantities with their discrete counterparts [9,10]: In two dimensions, by the Gauss-Bonnet theorem the Einstein-Hilbert action the Einstein-Hilbert action d 2 x g(x)R(x) is a topological invariant, which makes such a theory classically trivial since there are no equations of motion.Regge [2] gave a beautiful proof of this theorem in terms of the deficit angle.The sum over the deficit angles in two dimensions is proportional to the Euler characteristic, namely i δ i = 4π(1 − g).The corresponding term in the gravitational action can therefore be dropped.If one keeps the area A fixed to its initial value then classically dynamics can only arise from the R 2 -interaction term.Such a term was used in three-and higher-dimensional studies to cure the problem of unboundedness of the gravitational action [11].
For triangulated surfaces the Euler relation reads where N 0 , N 1 , and N 2 denote the number of sites, links and triangles, respectively.For triangulations without boundary we also know that a link is shared by two triangles, resulting in the relation From these two relations one can derive two more, namely 3, which will be useful later.
For each triangle there is a one-to-one correspondence between the square of the link lengths and the components of the metric.Denoting by g µν (i) the components of the metric tensor for the i th triangle, and by q i+µ,i+ν , q i,i+µ , and q i,i+ν the square of its three edge lengths, one can derive the relation g µν (i) = 1/2 [q i,i+µ + q i,i+ν − q i+µ,i+ν ].In classical Regge calculus one starts with the action principle and derives the equations of motion, one for each link.The link lengths have to be adjusted to satisfy those equations in order to be a classical solution.The connectivity of the edges, in simplicial topology called the incidence matrix, is usually fixed from the very beginning by the simplicial decomposition of the manifold under consideration.
In quantum Regge calculus the technical aspects are similar, although the philosophy is quite different.Here we want to evaluate the functional integral in equation ( 1) by, e.g., Monte Carlo (MC) simulation methods.In principle, the integral has to be extended over all metrics of all possible topologies, but commonly one restricts oneself to a specific topology, typically the sphere or the torus (the latter corresponding to periodic boundary conditions).The integral over the metric is replaced by an integral over the square of the link lengths.An important ingredient in the functional-integral method is the appropriate measure which is not even known in the continuum (this is a dramatic difference to path-integral quantization of particle mechanics).The most popular measure is DeWitt's supermetric [8], a distance functional on the space of metrics.It was used by Polyakov in his famous string solution [12].Because in two dimensions the measure is the primary source of the non-trivial dynamical content of the theory, its correct transcription might be the key point for a proper formulation.Nevertheless, if the discretised DeWitt measure is still a local one, then one might argue on the basis of universality that other local measures, in between some reasonable bounds close to the DeWitt measure, will do as well.In fact, most simulation studies reported in the literature employ a simplified scale-invariant, so-called "computer" measure where q ij = l 2 ij .The function F ({q ij }) takes on the value one if the update proposals for the link lengths do not violate the triangle inequalities, and it is zero otherwise.The parameter serves to suppress very thin triangles by generalising the triangle inequalities to a (still scale invariant) form (l 1 + l 2 ) (1 − )l 3 and (l 1 − l 2 ) (1 + )l 3 , which makes the algorithms somewhat faster.
Collecting the transcriptions from the continuum to the simplicial Regge approach, the lattice analogue to equation ( 1) is therefore given by where the abbreviation R 2 i = δ 2 i /A i was used.If "matter fields" are represented in a bold simplification by Ising spins σ i = ±1, their energy and coupling to the geometry is usually modelled by where the spins σ i are located at the vertices i of the lattice.Here the volume A ij associated with a link l ij is defined as Unfortunately the results of numerical investigations using Regge calculus have been quite disappointing so far.Using the commonly employed dl/l measure on a Regge lattice, no change in the phase transition of an Ising model coupled to gravity was observed [9,10], the critical exponents remained in the flat space Onsager universality class.Still, there is the hope that with a different measure or a different spin coupling to gravity one can reproduce the modified critical exponents predicted by the KPZ/DDK approach discussed in the next subsection [13,14].Measurements of the scaling properties of pure gravity, such as the string susceptibility exponent γ str , have themselves given rise to some disagreement in numerical investigations of the Regge calculus approach [15].

Dynamical triangulations and quadrangulations
An apparently more promising candidate for the construction of a consistent theory is the dynamical tessellations approach where all edge lengths of the simplicial building blocks are kept fixed and equal, but the connectivities are allowed to vary locally.For the Euclidean case in two dimensions, such an ensemble is commonly taken as the set of all gluings of equilateral triangles to a regular, usually closed surface of fixed topology, while counting each of the possible gluings with equal weight.Alternatively one may also consider quadrangulations as sketched in figure 1, where instead of triangles the simplicial building blocks are taken as quadrangles.The resulting random-surface model and its simplicial generalisation to higher dimensions are numerically tractable, for instance by Monte Carlo (MC) simulations.For two dimensions, the use of matrix models and generating-function techniques led to exact solutions for the cases of pure Euclidean gravity [16] and the coupling of certain kinds of matter, such as the Ising model [17,18], to the surfaces.Furthermore, the critical exponents governing phase transitions that the matter fields may exhibit are conjectured exactly from conformal field theory as functions of the exponents on regular lattices and the central charge C of the matter variables via the so-called KPZ/DDK formula [19] where ∆ is the original scaling weight and ∆ the "dressed" scaling weight upon coupling to gravity.The field-theory ansatz leading to equation ( 11) breaks down for central charges C > 1, an effect which has been termed the C = 1 "barrier", whereas discrete models of C > 1 matter coupled to dynamical triangulations or quadrangulations still appear to be well-defined.This mismatch of descriptions and its driving mechanism is still one of the less well understood aspects of the dynamical tessellations method.
For Monte Carlo simulations of 2D combinatorial dynamical triangulations or their dual φ 3 graphs, an ergodic set of updates is given by the so-called Pachner moves [20].An adaptation of these link-flip moves to simulations of quadrangulations proposed in [21] is shown in figure 2. Explicit counter-examples show, however, that these moves do not in general constitute an ergodic dynamics for simulations of dynamical quadrangulations.Introducing a second type of link- flip moves, a "two-link flip" (see figure 2), we recently constructed an algorithm for dynamical quadrangulations, which does not show any signs of ergodicity breaking [22][23][24].Analyses of autocorrelation times reveal, however, that the performance -as expected for a local algorithmis severely limited by slowing down near criticality.To alleviate this problem, we adapted the non-local "baby-universe surgery" method proposed in [25] for triangulations to quadrangulations [23,24].
For pure triangulations (no coupling to matter fields), independent realisations of this graph ensemble can also be generated more easily by a recursive insertion method proposed in [26].The dual graphs are planar, "fat" (i.e., orientable) φ 3 Feynman diagrams without tadpoles and self-energy insertions, which can be counted analytically by matrix model methods [3,16].

Exact solution for the Ising model on dynamical graphs
One remarkable result that emerged from studies of various statistical mechanical models coupled to two-dimensional quantum gravity is the exact solution of the Ising model in an external magnetic field [17].Even though this model appears quite complicated at first glance, the exact solution is more general than the one for the Ising model on regular lattices, for which the famous Onsager solution covers only the zero-field case.In discrete form the coupling of spin models to gravity may be interpreted from a statistical mechanics point of view as a special kind of annealed disorder in the form of random triangulations or quadrangulations, or their dual planar graphs.The partition function for the Ising model on a single graph G n with n vertices reads where σ i = ±1, and β = 1/T may be interpreted as inverse temperature and h as external magnetic field.Coupling to gravity means that this partition function is generalised by incorporating a sum over some class of graphs {G n }, Notice that the summations over the spin degrees of freedom ( {σ} ) and all graphs ( {G n } ) appear on the same footing.Viewed from the perspective of the Ising model, this represents annealed disorder for the spins.Generalisations to other "matter fields" are straightforward by replacing the Ising spins by, e.g., O(n) or Potts model spins or continuous field variables (with appropriate interaction terms).
The solution for the Ising model in [17] proceeded by first forming the grand-canonical partition function and then noting that this could be expressed as the free energy of a two-matrix model with N × N Hermitian matrices φ 1 and φ 2 .The coupling between the two types of fields is defined as c = exp(−2β).
The graphs of interest are generated as the Feynman diagrams of the "action" in equation ( 15), which is constructed in such a way that each edge or link of the graph carries the correct Boltzmann weight for Ising spins with nearest-neighbour interactions.Usually the N → ∞ limit is performed in order to pick out planar graphs (i.e., a well-defined spherical topology), but in a systematic 1/N expansion other topologies can be realised as well [18].In the other extreme N = 1, the model would generate all possible Feynman diagrams in a mean-field like manner.Since for general N > 1 the edges carry matrix indices, the graphs in question are sometimes called "fat" or ribbon graphs, while for N = 1 one speaks of "thin" or generic graphs.Due to the φ 4 terms in (15), the above matrix model generates four-valent, so-called φ 4 graphs.Alternatively, one could also consider (formal) φ 3 interactions and the resulting φ 3 graphs.
In the limit of large N the functional integral in (15) can be evaluated by saddle-point methods using the results of [27] to give where g(z) is defined by

Critical exponents
Equations ( 16) and ( 17) give an implicit solution for the grand-canonical partition function ( 14), from which the canonical Z n for any number of vertices n can be extracted by a series expansion.Note that, albeit only implicitly, this yields the exact answer also in an external magnetic field which is not available for the Ising model on regular lattices.There, Onsager's solution describes only the case h = 0, and with some extra effort the spontaneous magnetisation in zero field can be derived.
By analysing the solution ( 16), (17) in the thermodynamic limit n → ∞, Kazakov derived the exact critical exponents of the Ising model coupled to dynamical triangulations, in perfect agreement with the KPZ/DDK formula (11).To apply the latter, we need the relations between the standard critical exponents of the specific heat and magnetisation, C ∝ t and the correlation length exponent νd h = 3, where d h 2 is the fractal dimension of the random gravity graphs.Note that (i) α is negative, giving a third -order transition, and (ii) all exponents agree with the critical exponents of the three-dimensional spherical model.It is unclear whether the latter is purely coincidental or not.

Partition function zeros
Given the exact solution ( 16), (17), one can also try to understand the critical properties of the model by analysing the zeros of the canonical partition function.The idea that the zeros of the partition function could determine the phase structure of a spin model first appeared in Lee and Yang's work [28] who specifically considered zeros in the complex field plane -now commonly called Lee-Yang zeros.Somewhat later, Fisher [29] extended this idea also to zeros in the complex temperature plane -the so-called Fisher zeros.In both cases one studies how the non-analyticity characteristic of a phase transition appears from the partition function on finite lattices or graphs, which may be written as a polynomial for a lattice with m edges and n vertices, again with c = exp(−2β), y = exp(−2h).Lee-Yang and Fisher showed that the behaviour of the zeros of this polynomial in the complex y or c plane, in particular the limiting locus as m, n → ∞, determined the phase structure.Since then many applications and refinements of this approach have been reported in the literature [30][31][32].
For temperature-driven transitions, for simplicity in zero external field, the thermodynamic limit of the free energy on some class of lattices or graphs {G n } becomes where L is some set of curves, or possibly regions, in the complex c plane on which the zeros have support and ρ(c) is the density of the zeros there.
The general question of how to extract the locus of zeros analytically has been considered by various authors, notably Shrock and collaborators [31] for Ising and Potts models.It was first observed in [30] that such loci could be thought of as Stokes lines separating different regions of asymptotic behaviour of the partition function in the complex temperature or field planes.More recently, the case of models with first-order transitions has been investigated by Biskup et al. [33], who showed that the partition function of a statistical mechanical model defined in a periodic volume V and depending on some complex parameter such as c or y, can be written in terms of complex functions F l describing k different phases, where the various F l are the metastable free energies per unit volume of the phases, and the real part F l = F characterises the free energy when phase l is stable.The zeros of the partition function are then determined from the solutions of the equations These equations are thus in perfect agreement with the idea that the loci of zeros should be Stokes lines, since the zeros of Z lie on the complex phase coexistence curves F l = F m in the complex parameter plane.The specific Biskup et al. results apply to models with first-order transitions, but we are interested here in an Ising model with a third-order transition, so it might initially seem that these results were inapplicable.We are saved by the fact that when considered in the complex temperature plane the transition is continuous only at the physical (i.e., real) point itself (and possibly some other finite set of points).The determination of the locus of Fisher zeros for the Ising model on random graphs in the thermodynamic limit using the ideas of the previous section turns out to be rather straightforward, as we now describe.Since we wish to match F between the various solution branches to obtain the loci of Fisher zeros and F ∼ log(g(c)) for the Ising model on planar graphs, the equation which determines the loci of zeros in the thermodynamic limit is where the low-temperature solution g L (c) and the various high-temperature solutions g Hi (c) are given by solving g (z) = 0 in equation ( 17) for z.The resulting curve is shown in the c plane in figure 3 where it can be seen that in addition to the physical phase transition at c = exp(−2β) = 1/4, an unphysical transition with the same KPZ [19] exponents appears at c = −1/4.The interior of the curve is the ferromagnetic FM region and the exterior the paramagnetic PM and unphysical "O" phases, separated by cuts on the imaginary axis which we have not shown.
The diamonds plotted in figure 3 are generated from a series expansion of Z in equation ( 16), which is arrived at by reverting the expression for g(z) and substituting the resulting z(g) into equation (16).Earlier work reported in [34] obtained similar results at lower orders.The form of the expression for Z means that the contributions from each of the terms in equation ( 16) are proportional to each other [35], so the full series for Z can be generated from 1/2 log (z/g).
The loci of Fisher zeros are highly non-universal, and we also show the zeros on "thin", generic random φ 3 graphs for comparison in figure 4. Recall that these models can be thought of as the N → 1 limit of the matrix models, rather than the N → ∞ planar limit.Similar methods to those discussed above also serve in this case where one has mean-field behaviour [36].For the Ising model on thin graphs F ∼ log S, where S is the saddle point action for either the low-or high-temperature phase.The equivalent of equation ( 23) is then giving the locus plotted in figure 4. The locus of partition function zeros for Potts models and the locus of chromatic zeros are also accessible on thin graphs.
In summary, we have seen that an analytic determination of Fisher zeros for the Ising model on both fat and thin random graphs is possible, and that series expansions are easily obtainable.The general form of the solution also holds on (planar) random triangulations and φ 3 graphs and in non-zero field, so all of these can also be investigated.

Vertex models on quadrangulations
One of the most general classes of statistical mechanics models with discrete symmetry are 6-and 8-vertex models [37,38].For the 6-vertex model, the six possible arrow configurations and their Boltzmann weights ω i are sketched in figure 5.The partition function follows by summing the weights of all the allowed arrow configurations.Special cases can be mapped onto more well-known Ising and Potts models or graph colouring problems [38].For two-dimensional regular lattices, several of these vertex models have been solved exactly, yielding a very rich phase diagram with various transition lines as well as critical and multi-critical points [38].Hence, coupling this class of models to a fluctuating geometry of the dynamical triangulations type is of obvious interest, both as a prototypic model of statistical mechanics subject to annealed connectivity disorder and as a paradigmatic type of matter coupled to two-dimensional Euclidean quantum gravity.Recently, the use of matrix model methods similar to that sketched above for the Ising model led to a solution of the thermodynamic limit of a special 6-vertex model, the F model, coupled to planar φ 4 graphs [39].This model corresponds to a C = 1 conformal field theory, i.e., it lies on the boundary to the region C > 1, where the KPZ/DDK formula (11) breaks down.The locus of the F model is depicted in the phase diagram of figure 5 for a (static) square lattice where the model exhibits a Kosterlitz-Thouless (KT) phase transition at β c = ln 2 [37,38].The same type of transition is predicted on dynamical lattices, and in particular the critical coupling β c = ln 2 should agree with that on the square lattice [39].In addition, a special slice of the 8-vertex model could also be analysed via transformation to a matrix model [40].Due to intrinsic limitations, however, the analytical approach cannot reveal the behaviour of the matter-related observables and the details of the occurring phase transition or the fractal properties of the graphs such as, e.g., their fractal dimension d h .
Since vertex models are generically defined on four-valent lattices, instead of considering the more common dynamical triangulations or the dual planar, "fat" (i.e., orientable) φ 3 graphs, one has to work in numerical studies with the more intricate ensemble of dynamical quadrangulations or the dual φ 4 Feynman diagrams.To update the arrow configurations, in [41] the loop-cluster algorithm [42] was employed, slightly modified for the case of simulations on random lattices.Among the most easily measurable quantities are the internal energy U and the specific heat C v .The observed non-scaling of C v with system size (see figure 6) is a first evidence for the expected KT-like transition.By comparing the estimates of U = 0.333 355 (11) and C v = 0.2137 (12)  for the infinite square lattice [37,38] of U = 1/3 and C v = 28(ln 2) 2 /45 ≈ 0.2989, one is led to the conjecture that the critical internal energy of the F model is not affected by the coupling to random graphs.As can be seen in figure 6, this is specific to the critical point, where the curves for the two lattice types cross, and probably indicates an additional common symmetry at criticality.When the vertex model is coupled to quantum gravity, we expect a renormalisation or dressing of the critical exponents as prescribed by the KPZ/DDK formula (11), which should also marginally apply to the present limiting case C = 1.Assuming standard scaling relations the weights ∆ of the energy operator and ∆ P associated with the spontaneous staggered polarisation P 0 are related to the usual critical exponents by α where trivially d h = 2 for a regular lattice.For the expected infinite-order KT-like phase transition individual exponents are not properly defined but the exponent ratios entering finite-size scaling, still have a well-defined meaning.With C = 1, the KPZ/DDK formula (11) then predicts for the dressed exponents on dynamical random graphs (where d h ≈ 4) To check the prediction for γ/d h ν one can consider the FSS of the staggered polarisability χ (analogous to a susceptibility) at its maxima for finite graphs or at the transition point β c = ln 2. By analogy to the square-lattice case [43], one expects a FSS form including a leading effective correction term, For the square-lattice model one has ω χ = 2, whereas for the random-graph model the correction exponent is not known.Asymptotically both FSS sequences are expected to lead to the same exponents.Unfortunately, this is not at all obvious in the presence of large correction effects for the accessible graph sizes (recall the large fractal dimension d h ≈ 4 for dynamical lattices), and in particular analyses of the maxima data turned out to be very intricate [41].However, assuming a vanishing leading exponent and fitting χ(N 2 ) at β c to a purely logarithmic increase yields high-quality fits as demonstrated in figure 7 and thus verifies the prediction (26).
For the spontaneous polarisation P 0 (analogous to a magnetisation), the FSS ansatz can be taken similarly as leading at β c to the estimate β/d h ν = 0.469 (15), which is again consistent within error bars with the prediction (26).
The fractal graph properties can be characterised by several quantities.A particular useful one is the fraction of loops of length two which as a function of β exhibits a peak at β 0 = 0.6894 (54)  (cf.figure 7), in good agreement with β c = ln 2 ≈ 0.693.This observable, which clearly reflects the matter back-reaction on the graphs, turned out to be much more suitable for locating β c than the more traditional quantities such as the peak location of the polarisability [41].The string susceptibility exponent γ s is defined through Z(N 2 ) ∼ e µ0N2 N γs−3

2
. By decomposing the graphs into a self-similar tree of "baby universes", the distribution of so-called minBUs ("minimal neck baby universes") of size B can be used to determine This method, originally introduced for triangulations or φ 3 graphs [25,44], has been generalised to φ 4 graphs [23,41].Pure φ 4 graphs yield γ s = −1/2 in agreement with universality.For the F model with central charge C = 1, where the scaling form has again to be augmented with logarithmic corrections, the estimates [41] are compatible with γ s = 0 for β ln 2 (critical phase) and γ s = −1/2 for β > ln 2 (ordered phase), in agreement with the KPZ/DDK conjecture.For the fractal dimension d h , analytical work yields conflicting predictions (4.83 or ∞ as C → 1 [45]).By a FSS analysis of the (geometrical) two-point correlation function of the graphs and of their mean extent we obtained d h = 4, independent of β [41].

Potts models on quenched φ 3 gravity graphs
In the rest of the paper we now turn attention to the quenched situation, where the quantum gravity framework is merely used to generate random graphs with a specific connectivity or coordination number distribution.The paradigm for studies of the effect of quenched, random disorder on universal properties of critical phenomena are uncorrelated, randomly distributed couplings [46][47][48][49].This includes ferromagnetic random-bond models as well as the physically very different case of spin glasses, where competing interactions complement disorder with frustration [47,[50][51][52][53][54].For a continuous phase transition in the idealised pure system, the effect of random bonds has been convincingly shown by renormalization group analyses as well as numerical investigations to be able to induce a crossover to a new, disorder fixed point [48,[55][56][57][58][59].The question thus arises whether quenched connectivity disorder can also lead to a new disorder fixed point.Numerical simulation studies of spin models on quenched lattices of Voronoï-Delaunay type in two and three dimensions, however, suggested this not to be the case [60].
Starting from a distribution of points in the plane, a Voronoï cell in two dimensions is defined as the region of the plane which is closer to a given vertex than to any other.The three-valent vertices where these cells meet and the cell edges make up the Voronoï diagram.Accordingly, the structure geometrically dual to the Voronoï diagram is the Delaunay triangulation.For regularly placed vertices one recovers the Wigner-Seitz elementary cells of solid state physics.If the vertices are chosen at random, the resulting Voronoï-Delaunay graph is referred to as Poissonian random lattice since the vertices can be considered as realisation of a Poisson point process [61,62].
In what follows we shall focus on the resulting variation of co-ordination numbers q i of the triangulation resp.loop lengths of the dual graph, neglecting the fact of differing edge lengths.The  distribution of co-ordination numbers for dynamical triangulations is shown in figure 8, where for comparison the Voronoï-Delaunay case is included as well.Two snapshots of the resulting graph structures are depicted in figure 9. From the Euler relation ( 6), the average co-ordination number is a topological invariant for a fixed number N 2 of triangles in two dimensions, given for spherical topology by [3] The variance of co-ordination numbers is defined as µ 2 ≡ q 2 i − q i 2 .It turns out that the random variables q i in general are not independently distributed, but are reflecting a spatial correlation of the disorder degrees-of-freedom in addition to the trivial correlation induced by the constraint (29).For nearest-neighbour vertices these correlations are approximately described by the Aboav-Weaire law [61], q m(q) = (6 − a)q + b , where q m(q) is the number of edges of the neighbours of a q-sided cell, and a and b are some parameters [62,63].

Harris and Harris-Luck criteria
In a seminal paper, Harris [51] employed phenomenological scaling theory to argue that for uncorrelated disorder a crossover to a new universality class should not occur for systems with a specific-heat exponent α < 0. It is now widely believed that also the converse is true, i.e., a crossover does occur for systems with α > 0 [55,56,64].In the marginal case α = 0, realised, e.g., by the Ising model in two dimensions, the regular critical behaviour is merely modified by logarithmic corrections [48].Similarly, for systems exhibiting a first-order phase transition in the regular case, the introduction of quenched disorder coupling to the local energy density can weaken the transition to second (or even higher) order [54].While this scenario has been rigorously established for the case of two dimensions and an arbitrarily small amount of disorder [52,53,65], the situation for higher-dimensional systems is less clear.For a variety of systems in three dimensions, however, sufficiently strong disorder has been shown numerically [66][67][68] to be able to soften the transition to a continuous one.
The relevance of randomness coupling to the local energy density crucially depends on how fast fluctuations of the local transition temperature induced by fluctuations of the random variables in a correlation volume die out as the critical point is approached.For independent random variables, this decay occurs with an exponent of d/2 in d dimensions.The comparison of this power with the inverse correlation length exponent 1/ν leads to Harris' celebrated relevance criterion d/2 < 1/ν or, assuming hyper-scaling to be valid, α = 2 − νd > 0 = α c [51,69].
Spatial correlations of the disorder degrees of freedom lead to a modification of the fluctuations present in "typical" patches of the random system with respect to the behaviour expected from the central limit theorem for independent random variables, which is implicitly presupposed by Harris' arguments.Such correlations for a random-bond model have been considered occasionally [70][71][72][73] and altered relevance criteria have been proposed [70,74].Luck [74] has considered a class of irregular systems not covered by the random-bond paradigm, namely that of quasi-crystalline or aperiodic structures, and formulated a generalised relevance criterion.Although he did not consider the systems with connectivity disorder such as the random graph models to be considered here, his reasoning should also apply to these cases.Measuring distances between two graph vertices by the unique number of links traversed in the shortest path connecting them, we consider a spherical patch P of radius R on a triangulation, containing B(R) vertices.Then, the fluctuations of the average co-ordination number in P , around its expected value J 0 = q = 1/N 0 i q i = 6(1 − 4(N 2 + 4) −1 ) in general decay in the limit R → ∞ of large patches as defining the wandering exponent ω of the considered graph type.Near criticality, the fluctuation σ ξ (J) of the average co-ordination number in a correlation volume induces a local shift of the transition temperature proportional to where µ 2 ≡ q 2 i − q i 2 .For the regular critical behaviour to persist, these fluctuations should die out as the critical point t = 0 is approached.This is the case when ω does not exceed the threshold value where in the second equality hyper-scaling was assumed to be applicable.This means that quenched correlated disorder with ω > ω c (ν) may be a relevant perturbation and a new type of critical behaviour could occur.By recasting equation (33), this happens for a given random graph type for For uncorrelated disorder with ω = 1/2, α c = 0 and the Harris criterion is recovered.
In [75] the wandering exponent ω was numerically determined by sampling the fluctuations defined in equation (32) for random graphs of increasing size N 2 (cf.figure 10) and fitting the resulting exponents ω(N 2 ) to the finite-size scaling (FSS) ansatz where θ is an a priori unknown correction exponent.This yields [75] ω ∞ = 0.7473(98) (dynamical triangulations) , with A = −0.73(37) and θ = 0.264 (70), suggesting that ω = 3/4 in this case.The criterion (34) then implies a relevance threshold of α c = −2, i.e., that the connectivity disorder of quantum gravity graphs should alter the critical behaviour of all known standard models.
The result for Voronoï-Delaunay lattices turned out to be well consistent with ω = 1/2 which would result from correlations decaying with a power larger than d = 2 (see also [70]).A direct inspection of the correlation function of co-ordination numbers indicated an even exponential decay [75].Thus, the relevance criterion (34) reduces to the Harris criterion, i.e., Voronoï-Delaunay connectivity disorder should be a relevant perturbation for models with specific-heat exponent α > 0.

Analytical considerations
Given the fact that several spin models interacting with annealed connectivity disorder of gravity type can be solved exactly (or at least the critical exponents can be predicted from the KPZ/DDK formula), it is tempting to look for analytic solutions also in the quenched case [76].In this case the disorder average has to be performed at the level of the free energy, where C ij is the connectivity matrix of the graphs.Clearly, the non-linear operation "ln" in between the two summations makes a direct exact evaluation very difficult.This is the typical situation for quenched disordered systems and one may resort to the well-known replica trick which represents the logarithm in the following form: where in the last identity the order of taking the quenched disorder average and the limit n → 0 was formally interchanged.The merit of this procedure is that [Z n ] av takes the form of an annealed average, albeit now for a system with n replicas, Notice that the spins of all n replica are interacting among each other via the same connectivity matrix (which for each of the graphs is different).Similar to a gauge field this mediates interactions between the replica.For random-bond systems the resulting expression looks formally similar with C ij replaced by the random couplings J ij .If one assumes that the J ij are independent Gaussian variables, the summation over disorder can be explicitly performed and generates explicit interactions between the replica which are usually treated in perturbation theory combined with renormalization group analyses.In the present case, one may more simply argue that (39) represents an annealed system with n matter fields and hence a total central charge of C tot = nC.If one uses this formally in the KPZ/DDK formula (11) (of course, possible problems with the C = 1 barrier are ignored) and then formally performs the n → 0 limit (in which the C = 1 barrier problem is apparently cured . . .), one arrives at the following dressing formula for the conformal weights in the quenched case [76]:
All simulations were done close to the transition point with the Wolff single-cluster update algorithm.Primary observables are the energy and magnetisation which were stored in time-series files.Using reweighting techniques [80] it is straightforward to compute all relevant quantities in the finite-size scaling regime, e.g., the specific heat C = β 2 N [ e 2 − e 2 ] av and susceptibility χ = β N [ m 2 − |m| 2 ] av , but also mixed quantities such as the derivative d ln[ |m| ] av /dβ [80].As a function of N the finite-size scaling behaviour of these quantities is expected to be where C reg is a regular background term, α, ν, and γ are the usual critical exponents, d h = 4 is the fractal dimension of the (pure) φ 3 graphs, and the f i (x) are various FSS functions with x = (β − β c )N 1/νd h being the scaling variable.The correction terms indicated by [1 + • • • ] become unimportant for sufficiently large system sizes N .From least-squares fits one then obtains the critical exponents listed in table 1.
Looking at the results in table 1 it is clear that the exponent estimates are different from the exact values for regular 2D lattices, giving a clear indication that the connectivity disorder of planar random graphs is a relevant perturbation in the renormalization group sense, similar to the situation for random-bond disorder.Even more, for the 2D Ising model (q = 2) the values are unambiguously different, while for random-bond disorder only rather subtle logarithmic modifications are expected which are difficult to observe in numerical studies [81,82].Our estimates for q = 2 are not incompatible with both the quenched and annealed KPZ values at the level of accuracy we have achieved, but those for q = 4 definitely match none of the possible theoretical predictions.The 10-state model has clearly non-trivial exponents, thus unambiguously indicating the expected softening effect.Remarkably, the estimated q = 10 values are a good fit to the theoretical quenched q = 4 prediction; they are certainly incompatible with the q = 4 annealed values.It is also noteworthy that the q = 10 measurements (and also the q = 4 quenched theory predictions) violate a supposedly general bound [83] for quenched systems, 1/νd h < 1/2.Hyper-scaling implies that α/νd h should be negative if the bound holds, which also is in clear conflict with our directly measured value for q = 10.The numerical estimates of 1/νd h for q = 2 and q = 4, on the other hand, are consistent with the bound.Whether the failure of the q = 10 model to observe the bound is a consequence of the technical details of the averaging procedure as suggested in [84] or a result of the long-range correlations in the disorder is still unclear.
In this context it is worth mentioning a closely related study [85] of the Ising model on quenched random graphs which formally can be characterised by a central charge d = −5.In this notation [76], our case corresponds to d = 0.Even though the simulated d = −5 graphs were much smaller and the statistics poorer, in [85] very good agreement was obtained with the appropriate generalisation of the quenched prediction (40).
Let us conclude this section with a brief remark on the second type of random lattices.According to the Harris criterion, connectivity disorder from Poissonian random lattices should be relevant for the q = 3 Potts model with α = 1/3 > 0. The FSS analysis presented in [86] yields, however, a thermal scaling exponent in very good agreement with that for the regular lattice model.This is remarkable, since connectivity disorder couples to the local energy density, such that a relevant perturbation is expected to predominantly show up in the energy-related exponents.Whether the observed small, but significant difference of the magnetic exponents indicates the onset of a crossover to a new universality class or is merely an effect of neglected corrections to scaling, has to be checked by a more careful scaling analysis including corrections, possibly augmented by simulations for even larger lattices.Furthermore, models with larger values of the specific-heat exponent α, such as the q = 4 Potts model or the Baxter-Wu model [38], which both have an exponent α = 2/3, might be good candidates to check whether a change of critical behaviour can be induced at all by the connectivity disorder of Poissonian random lattices.

Summary
Annealed and quenched quantum gravity graphs provide a rich laboratory for analytical and numerical investigations of statistical physics systems.In the annealed case, the KPZ/DDK formula translates the critical behaviour on regular lattices to that when the same spin model is coupled to fluctuating graphs.In some situations it is even useful to turn this argument around; due to the possibility of matrix model solutions it is sometimes easier to find an exact solution for the random graph case and then translate this back to regular lattices.In the quenched case, the random graph ensembles are welcome prototypes for connectivity disorder with well-defined and in part exactly known properties.

Figure 1 .
Figure 1.A section of planar random quadrangulation (in bold) and the dual φ 4 graph (dashed).

Figure 8 .
Figure 8.Comparison of the co-ordination number distributions P (q) of Poissonian Delaunay triangulations and dynamical triangulations in the limit N2 → ∞.

Figure 9 .
Figure 9. Snapshots of random Poissonian Delaunay triangulations (left) and dynamical triangulations (right) of spherical topology with N2 = 5000 triangles.The Voronoï resp.φ 3 graphs considered numerically are the geometric duals of the shown structures.

Figure 10 .
Figure10.Numerical estimate of the scaling of the average fluctuation of co-ordination numbers of triangulations of volume N2 = 500 000 for the two considered ensembles and fits (bold lines) to the expected functional form(32).
at β c = ln 2 obtained on large φ 4 random graphs (N 2 = 65 536) with the analytical critical values Internal energy U and specific heat Cv of the F model coupled to dynamical quadrangulations.