Spontaneous Magnetization In the Plane

arXiv:cond-mat/0007303v1 [cond-mat.stat-mech] 19 Jul 2000 Spontaneous magnetisation in the plane Geoff Nicholls Department of Mathematics Auckland University Private Bag 92019, Auckland New Zealand [email protected] May 27, 1996 Revised July 20, 2000 The Arak process is a solvable stochastic process which generates coloured patterns in the plane. Patterns are made up of a variable number of random nonintersecting polygons. We show that the distribution of Arak process states is the Gibbs distribution of its states in thermodynamic equilibrium in the grand canonical ensemble. The sequence of Gibbs distributions form a new model parameterised by temperature. We prove that there is a phase transition in this model, for some non-zero temperature. We illustrate this conclusion with simulation results. We measure the critical exponents of this off-lattice model and find they are consistent with those of the Ising model in two dimensions. Figure Captions Figure 1 (A) A state χ of the Arak process (B) The discontinuity set γ of (A). Figure 2 (A) A set of lines ℓ intersecting D (B) an admissible graph drawn on the set ℓ (C) one of the two colourings of D with discontinuity set given by the graph in (B). Figure 3 Updates in the Markov Chain Monte Carlo. Dashed and solid edges are exchanged by the moves, which are reversible. (A) Interior vertex birth and death (B) move a vertex, and (C) recolour a region by swapping a pair of edges. In an extra move, not shown, a small triangle may be created or deleted. Further move types are used to update boundary structures. Figure 4 Binder parameter Ud (see text), regressed with cubic polynomials. Curves correspond to distinct box-side lengths d. The maximum likelihood fit, constrained to intersect at a point, is shown. Error bars in this and all other graphs are 1σ. Figure 5 The Binder parameter data of Figure 4 rescaled with Ising critical exponents. The regression is a cubic polynomial. χ243−4 = 38.5 for the fit is acceptable. Figure 6 The magnetisation m̄d (T ), regressed with cubic polynomials. Figure 7 The magnetisation data of Figure 6 rescaled with Ising critical exponents. The regression is a quartic polynomial. The value of the χ2 statistic shows that the fit is a poor one. Figure 8 A selection of states equilibrated in a box of side d = 12 at temperatures below and above the estimated critical temperature Tc ≃ 0.6665(5). 1 Introduction The Widom-Rowlinson model, with two species of discs and hard-core interactions between discs of unlike species, is sometimes referred to as the “continuum Ising model”. However there is another continuum model which might share the title. In 1982 Arak [1] presented a stochastic process in the plane with realisations of the kind shown in Figure 1A. States are composed of a variable number of coloured non-intersecting random polygons. Remarkably, the normalising constant is available as an explicit function of the area and boundary length of the region in which the process is realised. We present rigorous results and simulation based measurements related to critical phenomena in a two dimensional “continuum Ising model” derived from the Arak process. There are few rigorous results for continuum models of interacting extended two dimensional objects. Moreover, relatively few Monte Carlo simulation studies have been made, perhaps on account of the complexity of the simulation algorithms required. The Widom-Rowlinson model has a phase transition [2]. Its critical exponents have been measured and put it in the Ising universality class [3]. Critical phenomena are known to occur in a range of related models with q ≥ 2 species and certain soft-core interactions [4, 5]. Where critical exponents have been measured [6] the universality class seems to be the class of the corresponding q−species Potts model. For single-species models rigorous existence results for phase transitions have been given only in certain restricted models having area interactions [7, 8]. In the model we consider the interface between black and white regions summarises the state in the same way that Peierls’ contours parameterise an Ising system. The energy associated with a state is proportional to the length of the interface. In contrast to the Ising model, the vertices of the polygon forming the interface take positions in the plane continuum. At a temperature T = 1, the model we consider corresponds to the Arak process. For this value of the temperature the partition function equals the normalising constant of the corresponding Arak process. At smaller values of the temperature we are no longer dealing with an Arak process. We no longer have a closed form for the partition function. However the model remains well defined, and two phases coexist at temperatures bounded away from zero. Besides this result, which we prove, we estimate the critical exponents of the temperature-modified Arak model, using Markov chain simulation to generate realisations of the process. Values (obtained by “data-collapsing”) are consistent with the corresponding critical exponents of the Ising model. This is in accord with what we expect from the hypothesis of universality, since the ground state of the temperature-modified Arak model is two-fold degenerate, and states are two dimensional. Although there is no high temperature limit for polygonal models (a class of models including the Arak process) consistent polygonal models might play this role (this point is made in [9]). We give no rigorously determined upper bound on the critical temperature, although it is clear, from our simulations, that the consistent Arak process has a single phase. 1 2 The Arak process We now define the Arak process, following [10]. A state is a colouring map χ : D → J from each point in an open convex set D ⊂ ℜ2 , onto a set J of possible colours. See Figure 1A. We write ∂D for the set of points in the boundary of D. We consider the simplest case, J = {black, white}, of two colours. (n) Let XD be the class of all finite subsets x of D ∪ ∂D. Let XD (n ≥ 1) (n) be the set of point-sets x composed of n points, so that XD = ∪∞ n=0 XD , with (0) XD = {φ} the subset x = φ containing no points. Let dxi be the element of area in D and length on ∂D. A measure dν(x) is defined on XD by dν(x) = dx1 dx2 . . . dxn . This is the measure of an independent pair of Poisson point processes of unit intensity, on the boundary and interior. Let ΓD (x) be the set of all “polygon graphs” γ which can be drawn on the point-set x, ie the set of all graphs which can be drawn in D with edges nonintersecting straight lines, with the points in x as vertices. All interior vertices must have degree 2 (they are V -vertices), and all boundary vertices degree 1 (I-vertices). γ is composed of a number of separate polygons which may be chopped off by the boundary. See Figure 1B. The space of all allowed polygon graphs is the union over vertex sets x of the polygon graphs of x: [ ΓD ≡ ΓD (x). x∈XD We define a measure on ΓD by dλ(γ) = κ(γ) dν(x(γ)), n Y 1 Y sin(ψi ), κ(γ) = e <i,j> ij i=1 (1) (2) for a pattern γ with vertices at x(γ) = (x1 , x2 . . . xn ). ψi is the smaller angle at vertex i in D, and the smaller angle made with the boundary tangent at xi for vertices on ∂D. The product over < i, j > runs over vertex pairs i, j connected by an edge in γ, with eij = |xi − xj | the length of the edge between vertices i and j. A counting measure is taken on the graphs of a fixed point set. The significance of κ is sketched at the end of this section. Arak’s probability measure on ΓD is PD {dγ} = 1 exp(−2L(γ)) dλ(γ), ZD (3) with L(γ) the summed length of all edges in γ, and ZD a normalising constant. Remarkably, ZD has a simple closed form [1, 10] (ie the model is solvable), ZD = exp(L(∂D) + πA(D)), 2 where L(∂D) and A(D) are respectively the perimeter length and area of D. Certain expectation values have been calculated (see [10, 11]). Some examples are given in Table 1. A colouring map χ : D → J is a function assigning a colour, black or white, to each point in D. See Figure 1A. Let a colouring map χ be given and let Bχ be the set of points x ∈ D with a black point, ie some y ∈ D such that χ(y) = black, in every ǫ-neighbourhood. Let Wχ be similarly defined for white points. Let γ(χ) = Bχ ∩ Wχ denote the discontinuity set of this colouring. For each polygon graph γ we consider two colouring maps χ : D → J each having discontinuity set γ(χ) = γ. The two distinct colourings of a given polygon graph are assigned equal probability, so the probability measure for colour maps is just PD {dγ(χ)}/2. The probability measure (3) has a number of beautiful properties, besides solvability. Striking are consistency and the Markov property. Consider an open region S of D with ǫ-neighbourhood (S)ǫ ⊂ D; the probability measure for events in S, given full information about χ on (∂S)ǫ , is independent of any further information about the state in D \ S. That is the Markov property. Next, let S be an open convex region S ⊂ D and let γS ∈ ΓS denote the restriction of a state γ ∈ ΓD to S. The probability measure for events simulated in D from PD {dγ} but observed in the subset S is equal to PS {dγ}, in other words PS {dγ} = PD {dγS }. That is consistency. The Arak process shares these properties with a much larger family of probability measures called the consistent polygonal models. See [10] for the general picture. We will now explain in brief how κ(γ) arises, following [10] closely. Consider a number of straight lines drawn in the plane. Let li = (ρi , φi ) where ρi is the perpendicular distance from the line to an origin and φi is the angle the line makes to the x-axis. The parameter space of li is L = [0, ∞) × [0, π). Let LD be that subset of L consisting of all lines intersecting D. Let dl = dρ dφ be Lebesgue measure of LD . Let LnD be the set of all line sets ℓ = {l1 , l2 . . . ln } made up of n lines, each in LD . In this parameterisation LD = ∪n LnD is the set of all sets of lines in the plane intersecting D, and dν̃(ℓ) = dl1 dl2 . . . dln is the element of measure of a line process in D, corresponding to a Poisson point process of unit intensity in LD . Referring to Figure 2, we define an admissible graph on a line set ℓ to be a graph with edges coinciding with lines in ℓ, such that each line in ℓ contributes a single closed segment of non-zero length to the graph. All interior vertices are V vertices, all boundary vertices are I vertices. The set of all admissible graphs which can be drawn on some line set in LD is identical to ΓD . Let γ be some legal graph drawn on the line set ℓ. Define a measure dλ̃(γ) = dν̃(ℓ) in ΓD using the line process as our base measure, and taking counting measure over the legal graphs of a line set. We now have two parameterisations of the graph: from its line set ℓ, or from its vertex set x. The authors of [10] have shown that dλ̃(γ) = dλ(γ), ie, κ(γ) arises as the Jacobian of the transformation between x and ℓ. 3 3 Properties of a temperature modified Arak process We choose to modify the measure (3), and consequently loose solvability. Consider a system of non-overlapping polygonal chains of fluctuating number, length and vertex composition, confined to a planar region D. The chains may be attached in some places to the boundary of D. The state is described by a graph γ ∈ ΓD . Micro-states are associated with elements of volume dλ̃(γ) in ΓD , so that in the Gibbs ensemble edge segments are isotropic in orientation (a rather unnatural choice). However, the Gibbs distribution QD {dγ} of this system is just the Arak distribution above, modified by the addition of a temperature parameter, as we now show. The Gibbs distribution QD {dγ} has a density, g(γ) say, with respect to dλ̃(γ), the line measure. The Boltzmann entropy of the system is Z H[g] = g(γ) ln(g(γ)) dλ̃(γ) ΓD In the grand canonical ensemble, the energy and dimension of the system state fluctuate about fixed average values. We suppose that the state energy E(γ) is given by the total length of the chains, E(γ) = cL(γ), with c a positive constant. The dimension of the vertex position vector x is dim(x) = 2ni + nb with ni (nb ) the number of interior (boundary) vertices in γ. Maximising the entropy subject to constraints on the mean energy and mean dimension of the state, we obtain the distribution of systems of chains, QD {dγ} = 1 exp(−cL(γ)/T ) q −ne dλ(γ), ZD T where T and q are Lagrange multipliers, and ne is the number of edges in γ (ne = ni + nb /2). Under the change of scale xi → q̂xi , the measure transforms as dλ(γ) → q̂ ne dλ(γ). We therefore set q = 1 without loss of generality. Setting c = 2 we obtain a “temperature-modified” Arak process PT,D {dγ} = 1 exp(−2L(γ)/T ) dλ(γ). ZD T (4) The function 2L(γ)/T is a potential, (ie ZD T is finite), at least when 0 ≤ T ≤ 1, and, by Theorem 8.1 of [12], the temperature-modified measure keeps the spatial Markov property of the Arak measure. W Let µB D (T ) be the mean proportion of D coloured black (and µD (T ) white), µB D (T ) = ET,D {A(Bχ )/A(D)} . The magnetisation of a state m(χ) = |A(Bχ ) − A(Wχ )|/A(D) measures the colour asymmetry in that state. In our simulations (reported below) we see a qualitatively Ising-like temperature dependence in the mean 4 magnetisation. We prove, in an Appendix, that there is long range order (ie phase coexistence) in magnetisation, at all sufficiently small temperatures. We have translated Griffiths’ version [13] of Peierls’ proof of phase coexistence in the Ising lattice model to this continuum case. B|W Let µD (T ) be the expected proportion of D coloured black given that the boundary is white, that is B|W µD (T ) = ET,D {A(Bχ )/A(D) ∂D ∩ Bχ = ∅} . Theorem For the temperature modified Arak process in an open convex region D ⊂ ℜ2 there exists a temperature Tcold, 0 < Tcold < 1 and a constant a, a > 0, such that B|W µD (Tcold ) ≤ 12 − a independent of the area A(D) of the region. Surgailis [9] has shown that, for an open convex set S ⊂ D, the thermodynamic limit D ր ℜ2 of PT,D {dγS } exists, for a class of measures including PT,D {dγ}, for all temperatures below some small fixed positive value. With the theorem above, B|W µB|W (T ) = lim 2 µD (T ) Dրℜ exists and satisfies µB|W (Tcold ) < 1/2 − a for some a > 0. Hence, there is phase coexistence at all temperatures T < Tcold . In fact it follows from the result stated in the Appendix that   1 1 4 8 µB|W (T ) ≤ , (5) + + 4π 2 z 3 z 2 z where z = (1/(πT ) − 1). Sketching the function of T on the right hand side of Equation (5), we see that Tcold > 0.18, though this bound is not at all sharp. Simulation (see below) shows that the model has a phase transition with critical temperature very close to T = 2/3. The proof of the theorem is in two parts. We are after an upper bound on the expected area coloured black. The area of black in a state with white boundaries is not more than the summed area of the polygons it contains, and is maximised when they are not nested. This observation leads to a simplified bound on the expected area coloured black, Equation (10). This first result is obtained by an obvious translation of the Griffiths calculation into the terms of a continuum process. In that case the next step, bounding the number ways a polygon can be drawn on a lattice of fixed size, using a fixed number of links, is straightforward. In the continuum, the analogous problem is to bound the volume of the parameter space of a polygon of fixed length, where volume is measured using λ̃, the line-based measure. The main difficulty lies in the fact that there are unbounded, but integrable, functions in the measure which arise, for example, when an edge length goes to zero; these would be absent if there were no polygon closure constraint; as a consequence the closure constraint may not be relaxed as simply as it is in the Ising case. 5 4 Simulation Results The probability measure PT,D {dγ} may be sampled using the Metropolis-Hastings algorithm, and Markov chain Monte Carlo. In our simulations we take D to be a square box of side length d. Note that the number of vertices is not fixed. Since the dimension of a state depends on the number of vertices in it, the Markov chain must make jumps, corresponding to vertex addition and deletion, between states of unequal dimension. Simulation algorithms of this kind are widely used in physical chemistry [14, 15] and statistics [16, 17]. Although there exist vertex birth and death moves sufficient for ergodicity, we allow a number of other moves in order to reduce the correlation time of the chain. See Figure 3. At each update we generate a candidate state γ ′ , by selecting one of the moves, and applying it to a randomly selected part of the graph. The candidate state becomes the current state (ie it is accepted) with a probability given by the Metropolis-Hastings prescription. Otherwise it is rejected and the current state is not changed. In this way a reversible Markov chain is simulated. The chain is ergodic, with equilibrium measure PT,D {dγ}. Full details of our algorithm, including explicit detailed balance calculations for all the Markov chain updates, are given in [11]. The sampling algorithm is quite complex, but because the model is solvable at T = 1, it is possible to debug the code, by comparing a range of estimated expectations with predicted values. In Table 1 we present a selection of system statistics at T = 1. Quantities in brackets are one standard deviation in the place of the last quoted digit. These data show how the simulation code was tested in comparisons using analytically derived expectation values. Let fˆ equal the average of some statistic f (χ) over an output sequence of length N , let ρf (t) equal the normalised autocorrelation (or ACF) of f at lag t and, for M > 0, let PM τf = 1 + 2 t=1 ρ(t) estimate the normalised autocorrelation time of f (χ) in the output, so that the variance of fˆ is estimated by τf var(f )/N . We used Geyer’s initial monotone indicator [18] to determine M , the lag at which the ACF is truncated. The asymptotic variance σρ2 of the ACF as t → ∞ was estimated and used as a consistency check on each measurement: the estimated ACF should fall off to zero smoothly, and at large lag should stay within 2σρ bounds of zero. As usual we cannot show the Markov simulation process has converged, but it is at least stationary. Run parameters for the measurements at T < 1 are summarised in Table 2. Autocorrelations reported are for T = 0.66, near the critical temperature. We estimate the integrated autocorrelation time τm of the state magnetisation, along with its standard error [19] and present these alongside the total run length. Referring to Table 2, the autocorrelation time is fitted within standard error by τm ∝ d4.6 . Our Metropolis Hastings algorithm is a local update algorithm and this places practical limits on the size of the largest system we can explore. We now report our measurements of the mean magnetisation, m̄d (T ) = ET,d {m(χ)}, and the Binder parameter Ud (T ) = 1 − ET,d {m(χ)4 } . 3ET,d {m(χ)2 }2 6 f (γ) L(γ) ne (γ) ni (γ) 2 ET =1,d {f (γ)} πd 4d + 4πd 4πd2 d L̂ n̂e n̂i 0.5 1.571(6) 5.13(2) 3.13(2) 1 3.14(1) 16.46(7) 12.47(6) 2 6.27(2) 58.1(3) 50.1(2) 4 12.55(2) 216.7(4) 200.6(4) 8 25.14(2) 836.2(5) 804.2(5) χ25 1.1 4.0 5.1 Table 1: Listed are a selection of estimates made from output at T = 1. Here d is the box side, and for a state γ, L(γ) is the total edge length, ne (γ) equals the number of edges, and ni (γ) equals the number of interior vertices. Quantities in brackets are one standard deviation and in the place of the last quoted digits. Under the scaling hypothesis, the various curves Ud (T ) indexed by d all intersect at a single T -value, the critical temperature [20], T = Tc say. A Bayesian estimate T̂c may be given for the intersection point. Let Û denote the ordered set of independent U -measurements we made (43 in all), let TU denote the ordered set of T values at which measurements were made, and let Σ̂U denote the ordered set of estimated standard errors for the measurements in Û . These data are represented by the error bars in Figure 4. Each measurement is an independent measurement. For each d = 6, 8, 12, 16, we model the unknown true curve Ud (T ) using a cubic Ud∗ (T ) = U ∗ + (T − T ∗ ) 2 X p a(d) p x . p=0 The parameterisation constrains the regression in such a way that the four curves intersect at a point (T ∗ , U ∗ ). We simulate the joint posterior distribution of the random variables a(6) , a(8) , a(12) , a(16) , T ∗ , U ∗ |Û , TU , ΣU , conditioning the slope to be negative in the region containing the data, and conditioning the lines to intersect at a point, but otherwise taking an improper prior equal to a constant for all vectors of parameter values. Again MCMC simulation was used. The marginal posterior distribution of T ∗ is very nearly Gaussian. Our estimate of the critical temperature is then T̂c = 0.6665(5). The quoted standard error is the standard deviation of T ∗ in its marginal posterior distribution. 7 1 # Updates ×106 16 τ̂m ×106 0.00181(7) 2 40 0.018(1) 4 6400 0.41(1) 6 16000 2.4(3) 8 128000 9.2(2) 12 300000 59(3) 16 300000 240(40) d Table 2: Listed are run parameters for simulations at T = 0.66, a temperature close to the measured critical temperature. An update is a single pass through the Metropolis-Hastings propose/accept simulation sequence. Measurements made at the same d value, but different temperatures, are based on the same number of updates. The Bayesian inference scheme used to estimate Tc above is attractive for several reasons. Above all it quantifies the uncertainty in our estimate of Tc , taking full account of the complex constraints applying in the regression (though taking no account of possible errors due to violations of scaling). The sensitivity of the outcome to the orders of the regressing polynomials was explored. The chosen orders were the smallest that gave an acceptable likelihood. The posterior mode, which is the maximum likelihood estimate for Tc , on account of our flat prior, occurs at T ∗ = 0.6663. Metric factors weight the mass of probability in the posterior distribution only slightly away from the maximum of the likelihood. Because the energy has a discrete two-fold symmetry, and states are two dimensional, we expect the model to lie in the universality class of the Ising model. Finite size scaling under the scaling hypothesis leads to a system size dependence of the form [20] m̄d (τ ) = d−β/ν g(d1/ν τ ) Ud (τ ) = f (d1/ν τ ) with f and g unknown functions, τ the reduced temperature (T /Tc − 1), and β and ν critical exponents. If we plot Ud (τ ) or dβ/ν m̄d (τ ) against d1/ν τ , we expect to see no significant dependence on system size d for τ near zero. Using the Ising critical exponents ν = 1 and β = 0.125 and our estimate T̂c for the critical temperature, we show, in Figures 5 and 7, the maximum likelihood fit to the transformed data. The transformed Ud -data lies on a smooth curve. The transformed m̂d -data does not give a satisfactory χ2 (all of the misfit comes from points at T > Tc ), but this is to be expected: we are seeing scaling violations (a satisfactory fit to a quartic can be obtained (χ229−5 = 30) by dropping points at large T from the d = 6 and d = 8 data). If this is so, then the critical exponents of the Ising model the temperature dependent Arak process are equal at the 8 precision of our simulation analysis. Sample realisations from the model, taken at temperatures around the critical temperature are shown in Figure 8. Acknowledgements It is a pleasure to thank Bruce Calvert (Mathematics, Auckland University) for his advice and ideas. Appendix: long range order We now give the proof. Condition on a white boundary. There can be no boundary vertices. Let ΓW D be the subspace of ΓD of polygon graphs with no boundary vertices. Let ΘD be the subspace of ΓW D of graphs made up of just one polygon. Each point in ΘD corresponds to a single polygon, lying wholly in D. We begin by proving the inequality Equation (10) below. Among states built from a given set of polygons, with no edge connected to the boundary, the black area is largest when the polygons are arranged so that none are nested. It follows that the area of black in a state χ with a white boundary is less than or equal to the sum of the areas of all the polygons in that state. The area of a polygon θ of perimeter length L(θ) is smaller than the area of a circle with the same perimeter, so A(θ) < L(θ)2 /4π and X L(θ)2 A(Bχ ) ≤ . (6) 4π θ⊂γ(χ) We want to take expectations of either side of Equation (6) so we clear γ from the domain of the sum, using Z X f (θ) δ(θ ⊂ γ(χ)) dν(x(θ)). f (θ) ≡ θ⊂γ(χ) ΘD δ(θ ⊂ γ(χ)) puts a delta function at each point in Θ corresponding to a polygon in γ. Each of these is a product of delta functions in D for the vertices of θ to coincide with those of γ, with an indicator function for the edge connections to coincide. x(θ) is the set of vertex coordinate variables of the polygon θ. Now take the expectation of A(Bχ )/A(D) over patterns χ in ΓW D . We have Z L(θ)2 B|W µD ≤ E{ δ(θ ⊂ γ(χ)) | ∂D ∩ Bχ = ∅ } dν(x(θ)). ΘD 4πA(D) The expectation of the delta function is by definition R δ(θ ⊂ γ) × e−2L(γ)/T dλ(γ) . E{ δ(θ ⊂ γ) | ∂D ∩ Bχ = ∅ } = ΓW D R e−2L(γ)/T dλ(γ) ΓW D (7) Simplify the denominator by restricting the integral to those graphs to which the polygon θ could be added without intersecting an edge of a polygon already in place. That is, if ΓθW D ≡ {γ ∈ ΓW D : γ ⊃ θ} 9 is the set of polygon graphs containing the polygon θ, then [ {γ \ θ} Γ̃θW D ≡ γ∈ΓθW D is the sub-domain of interest. We have Z Z −2L(γ)/T e dλ(γ) ≥ e−2L(γ)/T dλ(γ) (8) Γ̃θW D ΓW D We now turn to the numerator of Equation (7). Carrying out the integration over vertices in θ using the δ-function, Z Z −2L(γ)/T δ(θ ⊂ γ) × e dλ(γ) = e−2L(γ)/T κ(θ) dλ(γ \ θ) ΓθW D ΓW D = κ(θ) e−2L(θ)/T Z e−2L(γ)/T dλ(γ), (9) Γ̃θW D since γ does not contain θ in the second line. Substituting with (8) and (9) in (7), and cancelling, E{ δ(θ ⊂ γ) | ∂D ∩ Bχ = ∅ } ≤ κ(θ) e−2L(θ)/T , and consequently, B|W µD Z 1 ≤ L(θ)2 e−2L(θ)/T dλ(θ). 4πA(D) ΘD In close analogy with Griffiths’ proof, we obtain Z  Z ∞ 1 B|W 2 −2b/T δ(b − L(θ)) dλ(θ) db b e µD ≤ 4πA(D) 0 ΘD (10) The integral over b is an integral over polygon perimeter lengths. The problem is now to bound the integral over ΘD without introducing more than one factor of A(D), or too rapidly growing a function of b. This is done by the following (n) Lemma. Let ΘD be the subset of ΘD of polygons with n vertices. Lemma Let Jn ≡ Z δ(b − L(θ)) dλ(θ) (n) (11) ΘD so that Z δ(b − L(θ)) dλ(θ) = ΘD ∞ X Jn n=3 in (10). Then Jn ≤ A(D)n2 (n − 1) and consequently Z (2π)n−1 bn−3 , (n − 2)! δ(b − L(θ)) dλ(θ) ≤ (2π)2 A(D)(4 + 2πb)2 e2πb . ΘD 10 (12) Proof of the Lemma: Start with Jn defined in Equation (11). Use a standard labelling with x1 the variable corresponding to the vertex in θ with the smallest x-coordinate, (smallest y-coordinate in case of ties) and vertex number increasing clockwise around θ. In the first step we break the polygon (n) at x1 to make a chain. Consider the set Θ̃D of distinct non-intersecting chains θ̃ of n edges linking n + 1 vertices, labeled with variables x1 to xn+1 . All the vertices in a chain lie entirely to the right of the first vertex (or directly above). (n) (n) Polygons are chains, ΘD ⊂ Θ̃D , since the first and last vertices in a chain may coincide. Transform variables from {xi }ni=1 to {x1 , {ei }ni=1 }, where ei is a Cartesian vector with origin xi corresponding to the edge from the i’th to the (i + 1)’th vertex. When we switch to integrating over chains, we constrain e1 + e2 + . . . + en to be zero, so that the polygon closes. Equation (11) becomes Z de1 de2 . . . den Jn ≤ δ(b − L(θ̃)) δ (2) (Σk ek ) dx1 , (n) e 1 e 2 . . . en Θ̃D with ei ≡ L(ei ) and using sin(ψi ) ≤ 1. The integrand is unbounded. We partition the space into regions, and impose the constraints b = L(θ̃) and e1 + e2 + . . . + en = 0 by integration over different variables in each region. For any particular region, the variables eliminated by the constraints are chosen so that the integrand is bounded in that region. Our second step then is to fix, by an integration in some dei , the closure constraint. We will need to be able to bound below the length of at least one edge of the chain. So define (n) Θ̃D,−ǫ (n,i) Θ̃D,ǫ (n) = { θ̃ ∈ Θ̃D = { θ̃ ∈ (n) Θ̃D | abs(L(θ̃) − b) ≥ ǫ } | abs(L(θ̃) − b) < ǫ, ei ≥ (b − ǫ)/n }, with i ∈ {1, 2, . . . , n}, and ǫ a small positive constant, 0 < ǫ < b, depending (n,i) on b. Each chain in Θ̃D,ǫ has the property that its ith edge has length at least (b − ǫ)/n. Any chain, with n edges and a total length differing from b by (n,i) not more than ǫ, must have such an edge. The sets Θ̃D,ǫ , i = 1, 2 . . . n are (n) (n) (n) not disjoint, but combine with Θ̃D,−ǫ to cover Θ̃D . Chains in ΘD,−ǫ will not contribute to the integral. It follows that Jn ≤ n Z X i=1 δ(b − L(θ̃)) δ (2) (Σk ek ) (n,i) Θ̃D,ǫ n ≤ n X (b − ǫ) i=1 Z δ(b − L(θ)) (n,i) ΘD,ǫ de1 de2 . . . den dx1 e 1 e 2 . . . en de1 de2 . . . de−i . . . den dx1 . e1 e2 . . . e−i . . . en (13) where a −i subscript indicates that element is left out of a product or sum. (n,i) ΘD,ǫ is the set of polygons with a long ith edge (that is, the set of chains in (n,i) Θ̃D,ǫ with xn+1 = x1 ). We have carried out the integral dei δ (2) (Σk ek ) and used the bound on ei . 11 The third step is to eliminate an edge length parameter, using b = L(θ̃), the length constraint. Let φi denote the angle made by edge ei to a fixed direction in the plane. In polar coordinates Equation (13) is n Z n X Jn ≤ δ(b − L(θ)) de1 dφ1 . . . de−i dφ−i . . . den dφn dx1 . (14) (b − ǫ) i=1 Θ(n,i) D,ǫ For the polygon to close ei sin(φi ) = − n X ek sin(φk ), (15) n X ek cos(φk ), (16) k=1 k6=i ei cos(φi ) = − k=1 k6=i and consequently b − L(θ) = b − n X ek − ek cos(φk − φi ). k=1 k6=i Integrating dej for some j may lead to an unbounded integrand. In order to (n,i) control this, we partition ΘD,ǫ on its angle variables. Let (n,i,j) ΘD,ǫ (n,i) = { θ ∈ ΘD,ǫ | π 3π < |φj − φi | < } 2 2 (n,i,j) A polygon in ΘD,ǫ has the property that the jth edge “turns back” from the direction of the long ith edge. There must be at least one such edge for the (n,i,j) polygon to close. The sets ΘD,ǫ , j = 1, 2 . . . n, j 6= i are not disjoint but their (n,i) union covers ΘD,ǫ . From Equation (14) n Jn ≤ n n XX (b − ǫ) i=1 j=1 j6=i Z δ(b − L(θ)) de1 dφ1 . . . de−i dφ−i . . . den dφn dx1 . (n,i,j) ΘD,ǫ We may now apply the integral dej to the delta-function δ(b − L(θ)). We transform from e, φ to e′ , φ′ where φ′k = φk and e′k = ek for 1 ≤ k ≤ n, k 6= j, and φ′j = φj and e′j = ej − ej cos(φj − φi ). The Jacobian of the full transformation e, φ → e′ , φ′ is just J −1 (e, φ → e′ , φ′ ) = ∂e′j ∂ej = 1 − cos(φj − φi ) − ej sin(φj − φi ) Repeated use of Equations (15) and (16) gives − sin(φj − φi ) ∂φi = , ∂ej ei 12 ∂φi . ∂ej (17) in Equation (17) and then using π/2 < |φj − φi | < 3π/2, we have J −1 > 1. The angle partition was needed to control this function. We can replace δ(b − L(θ)) dej by one, and restrict the integration domain to polygons of length b, ie set ǫ = 0. We obtain the simplified bound n n n XX Jn ≤ b i=1 j=1 j6=i Z de1 dφ1 . . . de−j dφj . . . de−i dφ−i . . . den dφn dx1 . (n,i,j) (18) ΘD,ǫ=0 (n,i,j) The last step is to bound the integral in Equation (18). Enlarge ΘD,ǫ=0 to allow each variable to range independently over its full domain, keeping only the bound on total edge length, L(θ) = b, and requiring x1 to remain in D. This will include polygons with crossing edges and allow the polygon to overlap the border of D. The integral dx1 gives a factor A(D). Each angle variable ranges over 0 to 2π contributing (2π)n−1 . The edge integrals are over the (n − 2)-dimensional tetrahedron e1 + e2 + . . . e−j + . . . e−i + . . . + en ≤ b − b/n of volume less than bn−2 /(n − 2)!. Combining these factors with a factor of (n − 1) from the sum over j, we obtain the bound on Jn given in the Lemma. This is the end of the proof of the Lemma. Equation (5) is obtained by evaluating the integral over b in Equation (10) with the bound from Equation (12), and the Theorem follows directly from Equation (5). 13 References [1] T. Arak. On Markovian random fields with a finite number of values. In 4th USSR-Japan Symposium on Probability Theory and Mathematical Statistics: Abstracts of Communications. Tiblisi, 1982. [2] D Ruelle. A phase transition in a continuous classical system. Phys. Rev. Lett., 27:1040–1041, 1971. [3] G Johnson, H Gould, J Machta, and LK Chayes. Monte Carlo study of the Widom-Rowlinson fluid using cluster methods. Phys. Rev. Lett., 79:2612– 2615, 1997. [4] JL Lebowitz and EK Lieb. Phase transitions in a continuum classical system with finite interactions. Phys. Lett., 39A:98–100, 1972. [5] HO Georgii and O Häggström. Phase transition in continuum Potts models. Commun. Math. Phys., 181:507–528, 1996. [6] R Sun, H Gould, J Machta, and LW Chayes. Cluster Monte Carlo study of multi-component fluids of the Stillinger-Helfand andWidom-Rowlinson type. Technical report, Physics Dept., Clark University, Worcester, USA, 2000. [7] JL Lebowitz, AE Mazel, and E Presutti. Rigorous proof of a liquid-vapor phase transition in a continuum particle system. Phys. Rev. Lett., 80:4701– 4704, 1998. [8] HO Georgii. Phase transition and percolation in Gibbsian particle models. In Conference Proceedings of Wuppertal, February 1999, Lecture Notes in Physics. Springer Verlag, 1999. [9] D Surgailis. The thermodynamic limit of Polygonal models. Acta Appl. Math, 22:77–102, 1991. [10] T. Arak, P. Clifford, and D. Surgailis. Point based Polygonal models for random graphs. Advances in Applied Probability, 25:348–372, 1993. [11] P Clifford and GK Nicholls. Simulating Polygonal shape models with data. Technical report, Oxford University Statistics Department, 1994. [12] T Arak and D Surgailis. Markov fields with polygonal realizations. Probability Theory Related Fields, 80:543–579, 1989. [13] RB Griffiths. Peierls’ proof of spontaneous magnetization of a two dimensional Ising ferromagnet. Phys. Rev. A, 136:437–439, 1964. [14] GE Norman and VS Filinov. Investigations of phase transitions by a Monte Carlo method. High Temperature, 7:216–222, 1969. Translation, Journal also known as High Temperature Research USSR. [15] D Frenkel. Advanced Monte Carlo techniques. In MP Allen and DJ Tildesley, editors, Computer simulation in chemical physics, volume C397 of Nato ASI Series. Kluwer Academic Publications, Dordrecht, 1993. [16] P. J. Green. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika, 82:711–732, 1995. 14 [17] C. J. Geyer and J. Møller. Simulation and likelihood inference for spatial point processes. Scandinavian Journal of Statistics, 21:359–373, 1994. [18] CJ Geyer. Practical Markov chain Monte Carlo. Statistical Science, 7:473– 511, 1992. [19] A Sokal. Monte Carlo methods in Statistical Mechanics. In Cours de Troisième Cycle de la Physique en Suisse Romande, Lausanne, 1989. [20] K Binder. Finite size scaling analysis of Ising model block distribution functions. Z. Phys. B, 43:119–140, 1981. B A Figure 1: (A) A state χ of the Arak process (B) The discontinuity set γ of (A). 15 A B C Figure 2: (A) A set of lines ℓ intersecting D (B) an admissible graph drawn on the set ℓ (C) one of the two colourings of D with discontinuity set given by the graph in (B). A B C Figure 3: Updates in the Markov Chain Monte Carlo. Dashed and solid edges are exchanged by the moves, which are reversible. (A) Interior vertex birth and death (B) move a vertex, and (C) recolour a region by swapping a pair of edges. In an extra move, not shown, a small triangle may be created or deleted. Further move types are used to update boundary structures. 16 0.7 d=6 d=8 d=12 d=16 T*,U* 0.6 0.5 U d 0.4 0.3 0.2 0.1 0 0.63 0.64 0.65 0.66 0.67 0.68 0.69 0.7 0.71 T Figure 4: Binder parameter Ud (see text), regressed with cubic polynomials. Curves correspond to distinct box-side lengths d. The maximum likelihood fit, constrained to intersect at a point, is shown. Error bars in this and all other graphs are 1σ. 17 0.8 0.7 0.6 Ud 0.5 0.4 0.3 0.2 0.1 0 −1 −0.8 −0.6 −0.4 −0.2 1/ν d 0 0.2 0.4 0.6 0.8 1 (T/Tc−1) Figure 5: The Binder parameter data of Figure 4 rescaled with Ising critical exponents. The regression is a cubic polynomial. χ243−4 = 38.5 for the fit is acceptable. 18 0.7 0.6 m d 0.5 0.4 0.3 0.2 0.1 0.63 0.64 0.65 0.66 0.67 0.68 0.69 0.7 0.71 T Figure 6: The magnetisation m̄d (T ), regressed with cubic polynomials. 19 1 0.9 0.8 0.7 0.5 d β/ν md 0.6 0.4 0.3 0.2 0.1 0 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 d1/ν(T/T −1) 0.4 0.6 0.8 1 c Figure 7: The magnetisation data of Figure 6 rescaled with Ising critical exponents. The regression is a quartic polynomial. The value of the χ2 statistic shows that the fit is a poor one. 20 T=0.62 0.64 0.66 0.68 0.70 0.72 Figure 8: A selection of states equilibrated in a box of side d = 12 at temperatures below and above the estimated critical temperature Tc ≃ 0.6665(5). 21