Spontaneous magnetisation in the plane

(1)

Spontaneous magnetisation in the plane

Geoff Nicholls

Department of Mathematics Auckland University Private Bag 92019, Auckland

New Zealand

[email protected] May 1996, Revised Sept 2000

The Arak process is a solvable stochastic process which generates coloured pat- terns in the plane. Patterns are made up of a variable number of random non- intersecting 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 pa- rameterised by temperature. We prove that there is a phase transition in this model, for some non-zero temperature. We illustrate this conclusion with sim- ulation 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.

arXiv:cond-mat/0007303 19 Jul 2000

Department of Mathematics, University of Auckland, Technical report #455

(2)

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 stud- ies 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 withq ≥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 correspondingq−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.

(3)

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 ⊂ <², 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.

Let X_D be the class of all finite subsets x ofD ∪∂D. For each n≥0, let X_D⁽ⁿ⁾ be the set of point-setsx={x₁, x₂, . . . xn}composed ofnpoints, so that X_D =∪^∞n=0X_D⁽ⁿ⁾, with X_D⁽⁰⁾ ={∅} the subset x=∅ of D ∪∂D containing no points. In the processes we define below, the number of pointsninxwill vary randomly from one realisation to another. Let dxⁱ be the element of area in D and length on ∂D. For each n ≥ 0, the element of volume at some point x∈X_D⁽ⁿ⁾ is equal todν(x) where

dν(x) =dx₁dx₂. . . dxⁿ

and dν(∅) = 1. Thus dν(x) is the measure, in X_D, 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 inD with edges non- intersecting straight lines, with the points inxas 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 xof the polygon graphs ofx:

Γ_D≡ [

x∈X_D

Γ_D(x).

We define a measure on Γ_D by

dλ(γ) = κ(γ)dν(x(γ)), (1)

κ(γ) = Y

<i,j>

1 e^ij

Yn

i=1

sin(ψⁱ), (2)

for a polygonal graphγ∈Γ_D with vertices at x(γ) = (x₁, x₂. . . xⁿ). In Equa- tion (1), ψⁱ is the smaller angle at vertex i for vertices in D, and the smaller angle made with the boundary tangent at xⁱ for vertices on ∂D. The prod- uct over < i, j > runs over vertex pairs i, j connected by an edge in γ, with e^ij = |xⁱ−x^j| the length of the edge between vertices i and j. A counting measure is taken on Γ_D(x). The significance ofκis sketched at the end of this section.

Arak’s probability measure on Γ_D is P_D{dγ}= 1

Z_D exp(−2L(γ))dλ(γ), (3)

withL(γ) the summed length of all edges inγ, andZD a normalising constant.

Remarkably,ZD has a simple closed form [1, 10] (iethe model issolvable), ZD = exp(L(∂D) +πA(D)),

(4)

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. LetW^χ 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 P_D{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 inS, 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 toS. The probability measure for events simulated in D from P_D{dγ} but observed in the subset S is equal to P^S{dγ}, in other words P^S{dγ} = P_D{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. Letli= (ρi, φi) whereρi is the perpendicular distance from the line to an origin and φi is the angle the line makes to thex-axis. The parameter space of li isL= [0,∞)×[0, π). LetLD

be that subset of L consisting of all lines intersecting D. Let dl = dρ dφ be Lebesgue measure of LD. Let LⁿD be the set of all line sets` ={l₁, l₂. . . ln} made up ofnlines, each inL^D. In this parameterisationLD=∪ⁿLⁿD is the set of all sets of lines in the plane intersectingD, and

d˜ν(`) =dl₁dl₂. . . dlⁿ

is the element of measure of a line process inD, corresponding to a Poisson point process of unit intensity inL^D. 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 areV vertices, all boundary vertices areI vertices.

The set of all admissible graphs which can be drawn on some line set inL_D 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 setx. The authors of [10] have shown thatdλ(γ) =˜ dλ(γ),ie, κ(γ) arises as the Jacobian of the transformation betweenxand`.

(5)

3 Properties of a temperature modified Arak process

We choose to modify the measure (3), and consequently loose solvability. Con- sider a system of non-overlapping polygonal chains of fluctuating number, length and vertex composition, confined to a planar regionD. The chains may be at- tached in some places to the boundary ofD. 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 Q_D{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 Q_D{dγ} has a density, g(γ) say, with respect to dλ(γ), the line measure. The Shannon entropy of the system is˜

H[g] = Z

ΓD

g(γ) ln(g(γ))d˜λ(γ)

In the grand canonical ensemble, the energy and dimension of the system state fluctuate about fixed average values. We suppose that the state energyE(γ) is given by the total length of the chains,E(γ) =cL(γ), withca positive constant.

The dimension of the vertex position vectorxis dim(x) = 2ni+nb withni (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,

Q_D{dγ}= 1

Z_DT exp(−cL(γ)/T)q⁻ⁿ^edλ(γ),

where T and q are Lagrange multipliers, and ne is the number of edges in γ (ne=ni+nb/2). Under the change of scalexi→qxˆ i, the measure transforms asdλ(γ)→qˆⁿ^edλ(γ). We therefore set q= 1 without loss of generality. Setting c= 2 we obtain a “temperature-modified” Arak process

P^T,_D{dγ}= 1

Z_DT exp(−2L(γ)/T)dλ(γ). (4) The function 2L(γ)/T is a potential, (ieZDT 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.

Letµ^BD(T) be the mean proportion ofDcoloured black (andµ^WD(T) white), µ^BD(T) =E^T,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

(6)

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.

Letµ^BD^|^W(T) be the expected proportion ofDcoloured black given that the boundary is white, that is

µ^BD^|^W(T) =E^T,D{A(B^χ)/A(D) ∂D∩B^χ=∅}.

Theorem For the temperature modified Arak process in an open convex region D ⊂ <² there exists a temperatureT_cold,0< T_cold<1 and a constanta,a >0, such that

µ^BD^|^W(T_cold)≤¹₂−a independent of the areaA(D)of the region.

Surgailis [9] has shown that, for an open convex set S ⊂ D, the thermo- dynamic limitD % <² of P^T,_D{dγ^S} exists, for a class of measures including P^T,_D{dγ}, for all temperatures below some small fixed positive value. With the theorem above,

µ^B^|^W(T) = lim

D%<²µ^BD^|^W(T)

exists and satisfiesµ^B^|^W(T_cold)<1/2−afor somea >0. Hence, there is phase coexistence at all temperaturesT < T_cold.

In fact it follows from the result stated in the Appendix that µ^B^|^W(T)≤ 1

4π² 1

z³+ 4 z²+8

z

, (5)

wherez= (1/(πT)−1). Sketching the function ofT on the right hand side of Equation (5), we see that T_cold >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 toT = 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.

(7)

4 Simulation Results

The probability measureP^T,_D{dγ}may be sampled using the Metropolis-Hastings algorithm, and Markov chain Monte Carlo. In our simulations we takeDto be a square box of side lengthd. 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 isrejectedand the current state is not changed. In this way a reversible Markov chain is simulated. The chain is ergodic, with equilibrium measure P^T,_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 atT = 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. The analytically derived expectation values given in the second line of the table come from [11]. They are derived using the particle representation of the Arak process given in [10]. Let ˆf equal the average of some statisticf(χ) over an output sequence of lengthN, letρ^f(t) equal the normalised autocorrelation (or ACF) off at lag tand, forM >0, letτf = 1 + 2PM

t=1ρ(t) estimate the normalised autocorrelation time off(χ) in the output, so that the variance of ˆf is estimated by τ^fvar(f)/N. We used Geyer’s initial monotone indicator [18] to determine M, the lag at which the ACF is truncated. The asymptotic variance σ²ρ 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 Ta- ble 2. Autocorrelations reported are forT = 0.66, near the critical temperature.

We estimate the integrated autocorrelation timeτ^mof 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 ∝d⁴^.⁶. 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) =

(8)

f(γ) L(γ) ne(γ) ni(γ) E^T=1,d{f(γ)} πd 4d+ 4πd² 4πd²

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)

χ²₅ 1.1 4.0 5.1

Table 1: Listed are a selection of estimates made from output atT = 1. Hered is the box side, and for a stateγ,L(γ) is the total edge length,ne(γ) equals the number of edges, andni(γ) equals the number of interior vertices. Quantities in brackets are one standard deviation and in the place of the last quoted digits.

E^T,d{m(χ)}, and the Binder parameter

U^d(T) = 1− E^T,d{m(χ)⁴} 3E^T,d{m(χ)²}².

Under the scaling hypothesis, the various curvesU^d(T) indexed bydall intersect at a single T-value, the critical temperature [20], T = T^c 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 TÛ denote the ordered set ofT values at which measurements were made, and let ˆΣÛ 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 curveUd(T) using a cubic

Ud^∗(T) =U^∗+ (T−T^∗) X2 p=0

a⁽p^d⁾T^p.

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⁽⁶⁾, a⁽⁸⁾, a⁽¹²⁾, a⁽¹⁶⁾, T^∗, U^∗|U, Tˆ ^U,Σ^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 ofT^∗ is very nearly

(9)

d # Updates ˆτ^m

×10⁶ ×10⁶

1 16 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)

Table 2: Listed are run parameters for simulations atT = 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 samed value, but different temperatures, are based on the same number of updates.

Gaussian. Our estimate of the critical temperature is then Tˆc = 0.6665(5).

The quoted standard error is the standard deviation ofT^∗in its marginal posterior distribution.

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 sensitiv- ity 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 forT^c, 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]

¯

md(τ) = d⁻^β/νg(d¹^/ντ) U^d(τ) = f(d¹^/ντ)

withf and g unknown functions, τ the reduced temperature (T/Tc−1), and β and ν critical exponents. If we plot Ud(τ) or d^β/νm¯d(τ) against d¹^/ντ, we expect to see no significant dependence on system sizedforτ near zero. Using the Ising critical exponentsν = 1 and β = 0.125 and our estimate ˆTc for the critical temperature, we show, in Figures 5 and 7, the maximum likelihood fit to the transformed data. The transformedU^d-data lies on a smooth curve. The

(10)

transformed ˆm^d-data does not give a satisfactoryχ²(all of the misfit comes from points atT > T^c), but this is to be expected: we are seeing scaling violations (a satisfactory fit to a quartic can be obtained (χ²₂₉₋₅= 30) by dropping points at largeT from thed= 6 andd= 8 data). If this is so, then the critical exponents of the Ising model the temperature dependent Arak process are equal at the 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

Here is the proof of the Theorem stated in Section 3. 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 Θ_Dcorresponds to a single polygon, lying wholly inD. 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 lengthL(θ) is smaller than the area of a circle with the same perimeter, soA(θ)< L(θ)²/4πand

A(B^χ)≤ X

θ⊂γ(χ)

L(θ)²

4π . (6)

We want to take expectations of either side of Equation (6) so we clearγ from the domain of the sum, using

X

θ⊂γ(χ)

f(θ)≡ Z

ΘD

f(θ)δ(θ⊂γ(χ))dν(x(θ)).

δ(θ⊂γ(χ)) puts a delta function at each point in Θ corresponding to a polygon inγ. Each of these is a product of delta functions in Dfor 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 ofA(Bχ)/A(D) over patternsχin ΓWD. We have µ^B_D^|^W ≤

Z

ΘD

L(θ)²

4πA(D) E{δ(θ⊂γ(χ))|∂D ∩Bχ=∅ }dν(x(θ)).

(11)

The expectation of the delta function is by definition E{δ(θ⊂γ)|∂D ∩Bχ=∅ }=

R

Γ^WDδ(θ⊂γ)×e⁻²^L⁽^γ⁾^/T dλ(γ) R

Γ^WDe⁻²^L⁽^γ⁾^/Tdλ(γ) . (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

Γ^θWD ≡ {γ∈Γ^W_D:γ⊃θ}

is the set of polygon graphs containing the polygonθ, then

˜Γ^θWD≡ [

γ∈Γ^θWD

{γ\θ}

is the sub-domain of interest. We have Z

Γ^WD

e⁻²^L⁽^γ⁾^/Tdλ(γ)≥ Z

˜Γ^θWD

e⁻²^L⁽^γ⁾^/Tdλ(γ) (8) We now turn to the numerator of Equation (7). Carrying out the integration over vertices inθ using theδ-function,

Z

Γ^WD

δ(θ⊂γ)×e⁻²^L⁽^γ⁾^/Tdλ(γ) = Z

Γ^θWD

e⁻²^L⁽^γ⁾^/Tκ(θ)dλ(γ\θ)

= κ(θ)e⁻²^L⁽^θ⁾^/T Z

˜Γ^θWD

e⁻²^L⁽^γ⁾^/Tdλ(γ), (9) sinceγdoes not containθ in the second line. Substituting with (8) and (9) in (7), and cancelling,

E{δ(θ⊂γ)|∂D ∩B^χ=∅ } ≤κ(θ)e⁻²^L⁽^θ⁾^/T, and consequently,

µ^BD^|^W ≤ 1 4πA(D)

Z

ΘD

L(θ)²e⁻²^L⁽^θ⁾^/Tdλ(θ).

In close analogy with Griffiths’ proof, we obtain µ^BD^|^W ≤ 1

4πA(D) Z _∞

0 b²e⁻²^b/T Z

ΘD

δ(b−L(θ))dλ(θ)

db (10)

The integral overbis an integral over polygon perimeter lengths. The problem is now to bound the integral over Θ_D without introducing more than one factor ofA(D), or too rapidly growing a function of b. This is done by the following Lemma. Let Θ⁽_Dⁿ⁾ be the subset of Θ_D of polygons with nvertices.

Lemma Let

Jⁿ ≡ Z

Θ⁽_Dⁿ⁾δ(b−L(θ))dλ(θ) (11)

(12)

so that Z

ΘD

δ(b−L(θ))dλ(θ) =X^∞

n=3

Jn

in (10). Then

Jn≤A(D)n²(n−1)(2π)ⁿ⁻¹bⁿ⁻³ (n−2)! , and consequently

Z

ΘD

δ(b−L(θ))dλ(θ) ≤ (2π)²A(D)(4 + 2πb)²e²^πb. (12)

Proof of the Lemma: Start with Jn defined in Equation (11). Use a standard labelling withx₁ 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 atx₁ to make a chain. Consider the set ˜Θ⁽_Dⁿ⁾of distinct non-intersecting chains θ˜ofn edges linking n+ 1 vertices, labeled with variables x₁ to xⁿ₊₁. All the vertices in a chain lie entirely to the right of the first vertex (or directly above).

Polygons are chains, Θ⁽_Dⁿ⁾ ⊂ Θ˜⁽_Dⁿ⁾, since the first and last vertices in a chain may coincide. Transform variables from {xⁱ}ⁿi=1 to {x1,{ei}ⁿi=1}, where ei is a Cartesian vector with origin xⁱ corresponding to the edge from the i’th to the (i+ 1)’th vertex. When we switch to integrating over chains, we constrain e₁+e₂+. . .+en to be zero, so that the polygon closes. Equation (11) becomes

Jn≤ Z

Θ˜⁽_Dⁿ⁾δ(b−L(˜θ))δ⁽²⁾(Σkek) de₁de₂. . . den

e₁e₂. . . en dx₁, withei≡L(ei) and using sin(ψi)≤1.

The integrand is unbounded. We partition the space into regions, and impose the constraintsb=L(˜θ) ande₁+e₂+. . .+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

Θ˜⁽_Dⁿ,⁾− = {θ˜∈Θ˜⁽_Dⁿ⁾ | abs(L(˜θ)−b)≥}

Θ˜⁽_D^n,i,⁾ = {θ˜∈Θ˜⁽_Dⁿ⁾ | abs(L(˜θ)−b)< , eⁱ≥(b−)/n},

with i ∈ {1,2, . . . , n}, and a small positive constant, 0 < < b, depending on b. Each chain in ˜Θ⁽_D^n,i,⁾ has the property that its ith edge has length at least (b−)/n. Any chain, with nedges and a total length differing frombby not more than , must have such an edge. The sets ˜Θ⁽_D^n,i,⁾, i = 1,2. . . n are not disjoint, but combine with ˜Θ⁽_Dⁿ,⁾− to cover ˜Θ⁽_Dⁿ⁾. Chains in Θ⁽_Dⁿ,⁾− will not

(13)

contribute to the integral. It follows that Jn ≤

Xn

i=1

Z

Θ˜⁽_D^n,i,⁾

δ(b−L(˜θ))δ⁽²⁾(Σkek) de₁de₂. . . den

e₁e₂. . . en dx₁

≤ n

(b−) Xn

i=1

Z

Θ⁽_D^n,i,⁾

δ(b−L(θ))de₁de₂. . . de₋i. . . den

e₁e₂. . . e₋ⁱ. . . eⁿ dx₁. (13) where a −i subscript indicates that element is left out of a product or sum.

Θ⁽_D^n,i,⁾ is the set of polygons with a long ith edge (that is, the set of chains in Θ˜⁽_D^n,i,⁾ with xn+1 = x₁). We have carried out the integral deiδ⁽²⁾(Σkek) and used the bound onei.

The third step is to eliminate an edge length parameter, usingb=L(˜θ), the length constraint. Letφⁱ denote the angle made by edgeei to a fixed direction in the plane. In polar coordinates Equation (13) is

Jⁿ ≤ n (b−)

Xn

i=1

Z

Θ⁽_D^n,i,⁾

δ(b−L(θ))de₁dφ₁. . . de₋ⁱdφ₋ⁱ. . . deⁿdφⁿdx₁. (14) For the polygon to close

eⁱsin(φⁱ) = − Xn

k=1 k6=i

e^ksin(φ^k), (15)

eⁱcos(φⁱ) = − Xn

k=1 k6=i

e^kcos(φ^k), (16)

and consequently

b−L(θ) =b− Xn

k=1 k6=i

ek−ekcos(φk−φi).

Integrating dej for some j may lead to an unbounded integrand. In order to control this, we partition Θ⁽_D^n,i,⁾on its angle variables. Let

Θ⁽_D^n,i,j, ⁾={θ∈Θ⁽_D^n,i,⁾ | π

2 <|φ^j−φi|<3π 2 }

A polygon in Θ⁽_D^n,i,j, ⁾has the property that thejth edge “turns back” from the direction of the long ith edge. There must be at least one such edge for the polygon to close. The sets Θ⁽_D^n,i,j, ⁾,j= 1,2. . . n, j6=iare not disjoint but their union covers Θ⁽_D^n,i,⁾. From Equation (14)

Jⁿ ≤ n (b−)

Xn

i=1

Xn

j=1 j6=i

Z

Θ⁽_D^n,i,j, ⁾

δ(b−L(θ))de₁dφ₁. . . de₋ⁱdφ₋ⁱ. . . deⁿdφⁿdx₁.

(14)

We may now apply the integral de^j to the delta-functionδ(b−L(θ)). We transform frome, φto e⁰, φ⁰ where φ⁰k =φ^k ande⁰k =e^k for 1≤k ≤n, k6=j, andφ⁰j =φ^j and

e⁰j=e^j−e^jcos(φ^j−φⁱ).

The Jacobian of the full transformatione, φ→e⁰, φ⁰ is just J⁻¹(e, φ→e⁰, φ⁰) = ∂e⁰j

∂ej

= 1−cos(φ^j−φⁱ)−e^jsin(φ^j−φⁱ)∂φⁱ

∂ej. (17) Repeated use of Equations (15) and (16) gives

∂φⁱ

∂e^j = −sin(φ^j−φⁱ)

eⁱ ,

in Equation (17) and then using π/2 < |φ^j −φⁱ| < 3π/2, we have J⁻¹ >

1. The angle partition was needed to control this function. We can replace δ(b−L(θ))de^jby one, and restrict the integration domain to polygons of length b,ieset = 0. We obtain the simplified bound

Jⁿ≤ n b

Xn

i=1

Xn

j=1 j6=i

Z

Θ⁽_D^n,i,j,=0⁾

de₁dφ₁. . . de₋^jdφ^j. . . de₋ⁱdφ₋ⁱ. . . deⁿdφⁿdx₁. (18)

The last step is to bound the integral in Equation (18). Enlarge Θ⁽_D^n,i,j,=0⁾ to allow each variable to range independently over its full domain, keeping only the bound on total edge length,L(θ) =b, and requiringx₁to remain inD. This will include polygons with crossing edges and allow the polygon to overlap the border ofD. The integraldx₁ gives a factorA(D). Each angle variable ranges over 0 to 2πcontributing (2π)ⁿ⁻¹. The edge integrals are over the (n−2)-dimensional tetrahedron

e₁+e₂+. . . e₋j+. . . e₋i+. . .+en≤b−b/n

of volume less than bⁿ⁻²/(n−2)!. Combining these factors with a factor of (n−1) from the sum overj, 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 overb in Equation (10) with the bound from Equation (12), and the Theorem follows directly from Equation (5).

(15)

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¨aggstr¨om. 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.

InConference 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. Proba- bility 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 Tildes- ley, editors,Computer simulation in chemical physics, volume C397 ofNato 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.

(16)

[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`eme 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.

A B

Figure 1: (A) A stateχof the Arak process (B) The discontinuity setγof (A).

(17)

B

A C

Figure 2: (A) A set of lines`intersecting D(B) an admissible graph drawn on the set` (C) one of the two colourings ofDwith 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.

(18)

0.630 0.64 0.65 0.66 0.67 0.68 0.69 0.7 0.71 0.1

0.2 0.3 0.4 0.5 0.6 0.7

T

U

d

d=6 d=8 d=12 d=16 T^*,U^*

Figure 4: Binder parameter Ud (see text), regressed with cubic polynomials.

Curves correspond to distinct box-side lengthsd. The maximum likelihood fit, constrained to intersect at a point, is shown. Error bars in this and all other graphs are 1σ.

(19)

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

d

^1/ν

(T/T

c

-1)

U

d

Figure 5: The Binder parameter data of Figure 4 rescaled with Ising critical exponents. The regression is a cubic polynomial. χ²₄₃₋₄ = 38.5 for the fit is acceptable.

(20)

0.63 0.64 0.65 0.66 0.67 0.68 0.69 0.7 0.71 0.1

0.2 0.3 0.4 0.5 0.6 0.7

T

m

d

Figure 6: The magnetisation ¯m^d(T), regressed with cubic polynomials.

(21)

-10 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 0.1

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

d

^1/ν

(T/T

c

-1)

d

β/ν

m

d

Figure 7: The magnetisation data of Figure 6 rescaled with Ising critical exponents. The regression is a quartic polynomial. The value of theχ² statistic shows that the fit is a poor one.

(22)

T=0.62 0.64

0.70

0.66 0.72 0.68

Figure 8: A selection of states equilibrated in a box of sided= 12 at temperatures below and above the estimated critical temperatureT^c'0.6665(5).