### 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

### 1 Introduction

The Widom-Rowlinson model, with two species of discs and hard-core interac- tions 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 realisa- tions 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 sim- ulation 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 mod-
els 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.

### 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.

Let *X** _{D}* be the class of all finite subsets

*x*of

*D ∪∂D. For each*

*n≥*0, let

*X*

_{D}^{(}

^{n}^{)}be the set of point-sets

*x*=

*{x*

_{1}

*, x*

_{2}

*, . . . x*

*n*

*}*composed of

*n*points, so that

*X*

*=*

_{D}*∪*

^{∞}*n*=0

*X*

_{D}^{(}

^{n}^{)}, with

*X*

_{D}^{(0)}=

*{∅}*the subset

*x*=

*∅*of

*D ∪∂D*containing no points. In the processes we define below, the number of points

*n*in

*x*will vary randomly from one realisation to another. Let

*dx*

*be the element of area in*

^{i}*D*and length on

*∂D. For each*

*n*

*≥*0, the element of volume at some point

*x∈X*

_{D}^{(}

^{n}^{)}is equal to

*dν(x) where*

*dν(x) =dx*_{1}*dx*_{2}*. . . dx*^{n}

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 in

*D*with edges non- intersecting 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}*≡* [

*x**∈**X*_{D}

Γ* _{D}*(x).

We define a measure on Γ* _{D}* by

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

*κ(γ) =* Y

*<i,j>*

1
*e*^{ij}

Y*n*

*i*=1

sin(ψ* ^{i}*), (2)

for a polygonal graph*γ∈*Γ* _{D}* with vertices at

*x(γ) = (x*

_{1}

*, x*

_{2}

*. . . x*

*). In Equa- tion (1),*

^{n}*ψ*

*is the smaller angle at vertex*

^{i}*i*for vertices in

*D, and the smaller*angle made with the boundary tangent at

*x*

*for vertices on*

^{i}*∂D. The prod-*uct over

*< i, j >*runs over vertex pairs

*i, j*connected by an edge in

*γ, with*

*e*

*=*

^{ij}*|x*

^{i}*−x*

^{j}*|*the length of the edge between vertices

*i*and

*j. A counting*measure is taken on Γ

*(x). The significance of*

_{D}*κ*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)

with*L(γ) the summed length of all edges inγ, andZ**D* a normalising constant.

Remarkably,*Z**D* has a simple closed form [1, 10] (iethe model is*solvable),*
*Z**D* = exp(L(∂D) +*πA(D)),*

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 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

*γ∈*Γ

*to*

_{D}*S. 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. Let*l**i*= (ρ*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* *l**i* is*L*= [0,*∞)×*[0, π). Let*L**D*

be that subset of *L* consisting of all lines intersecting *D. Let* *dl* = *dρ dφ* be
Lebesgue measure of *L**D*. Let *L*^{n}*D* be the set of all line sets*`* =*{l*_{1}*, l*_{2}*. . . l**n**}*
made up of*n*lines, each in*L** ^{D}*. In this parameterisation

*L*

*D*=

*∪*

^{n}*L*

^{n}*D*is the set of all sets of lines in the plane intersecting

*D, and*

*d˜ν(`) =dl*_{1}*dl*_{2}*. . . dl*^{n}

is the element of measure of a line process in*D, corresponding to a Poisson point*
process of unit intensity in*L** ^{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 are

*V*vertices, all boundary vertices are

*I*vertices.

The set of all admissible graphs which can be drawn on some line set in*L** _{D}* is
identical to Γ

*. Let*

_{D}*γ*be some legal graph drawn on the line set

*`. Define a*measure

*dλ(γ) =*˜

*d˜ν(`) in Γ*

*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*

_{D}*`, or from its vertex setx. The*authors of [10] have shown that

*dλ(γ) =*˜

*dλ(γ),ie,*

*κ(γ) arises as the Jacobian*of the transformation between

*x*and

*`.*

### 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 region*D. The chains may be at-*
tached 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 Γ*˜

*, so that in the Gibbs ensemble edge segments are isotropic in orientation (a rather unnatural choice). However, the Gibbs distribution Q*

_{D}

_{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 energy*E(γ) is*
given by the total length of the chains,*E(γ) =cL(γ), withc*a positive constant.

The dimension of the vertex position vector*x*is dim(x) = 2n*i*+*n**b* with*n**i* (n*b*)
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*_{D}*T* exp(−cL(γ)/T)*q*^{−}^{n}^{e}*dλ(γ),*

where *T* and *q* are Lagrange multipliers, and *n**e* is the number of edges in *γ*
(n*e*=*n**i*+*n**b**/2). Under the change of scalex**i**→qx*ˆ *i*, the measure transforms
as*dλ(γ)→q*ˆ^{n}^{e}*dλ(γ). 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*_{D}*T* exp(−2L(γ)/T)*dλ(γ).* (4)
The function 2L(γ)/T is a potential, (ie*Z**D**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.

Let*µ*^{B}*D*(T) be the mean proportion of*D*coloured black (and*µ*^{W}*D*(T) white),
*µ*^{B}*D*(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

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*µ*^{B}*D*^{|}* ^{W}*(T) be the expected proportion of

*D*coloured black given that the boundary is white, that is

*µ*^{B}*D*^{|}* ^{W}*(T) =E

^{T,}*D*

*{A(B*

*)/A(D)*

^{χ}*∂D∩B*

*=*

^{χ}*∅}.*

Theorem *For the temperature modified Arak process in an open convex region*
*D ⊂ <*^{2} *there exists a temperatureT*_{cold}*,*0*< T*_{cold}*<*1 *and a constanta,a >*0,
*such that*

*µ*^{B}*D*^{|}* ^{W}*(T

_{cold})

*≤*

^{1}

_{2}

*−a*

*independent of the areaA(D)of the region.*

Surgailis [9] has shown that, for an open convex set *S* *⊂ D, the thermo-*
dynamic limit*D % <*^{2} 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%<*^{2}*µ*^{B}*D*^{|}* ^{W}*(T)

exists and satisfies*µ*^{B}^{|}* ^{W}*(T

_{cold})

*<*1/2

*−a*for some

*a >*0. Hence, there is phase coexistence at all temperatures

*T < T*

_{cold}.

In fact it follows from the result stated in the Appendix that
*µ*^{B}^{|}* ^{W}*(T)

*≤*1

4π^{2}
1

*z*^{3}+ 4
*z*^{2}+8

*z*

*,* (5)

where*z*= (1/(πT)*−*1). Sketching the function of*T* 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 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.

### 4 Simulation Results

The probability measure*P*^{T,}_{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 ver-
tex 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 *γ** ^{0}*, 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

*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
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. 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
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

*τ*

*f*= 1 + 2P

*M*

*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 Ta-
ble 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}*∝d*

^{4}

^{.}^{6}. Our Metropolis Hastings algorithm is a local update algo- rithm 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) =

*f*(γ) *L(γ)* *n**e*(γ) *n**i*(γ)
E* ^{T}*=1

*,d*

*{f*(γ)}

*πd*4d+ 4πd

^{2}4πd

^{2}

*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)

*χ*^{2}_{5} 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,n**e*(γ) equals the
number of edges, and*n**i*(γ) 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(χ)*

^{4}

*}*3E

^{T,d}*{m(χ)*

^{2}

*}*

^{2}

*.*

Under the scaling hypothesis, the various curves*U** ^{d}*(T) indexed by

*d*all intersect at a single

*T*-value, the critical temperature [20],

*T*=

*T*

*say. A Bayesian estimate ˆ*

^{c}*T*

*may be given for the intersection point. Let ˆ*

^{c}*U*denote the ordered set of independent

*U*-measurements we made (43 in all), let

*T*

*denote the ordered set of*

^{U}*T*values at which measurements were made, and let ˆΣ

*denote the ordered set of estimated standard errors for the measurements in ˆ*

^{U}*U*. 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

*U*

*d*(T) using a cubic

*U**d** ^{∗}*(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*^{(6)}*, a*^{(8)}*, a*^{(12)}*, a*^{(16)}*, 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 of*T** ^{∗}* is very nearly

*d* # Updates ˆ*τ*^{m}

*×10*^{6} *×10*^{6}

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 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.

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 pos-
terior distribution.

The Bayesian inference scheme used to estimate *T**c* above is attractive for
several reasons. Above all it quantifies the uncertainty in our estimate of *T**c*,
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 for*T** ^{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]

¯

*m**d*(τ) = *d*^{−}^{β/ν}*g(d*^{1}^{/ν}*τ)*
*U** ^{d}*(τ) =

*f*(d

^{1}

^{/ν}*τ)*

with*f* and *g* unknown functions, *τ* the reduced temperature (T/T*c**−*1), and
*β* and *ν* critical exponents. If we plot *U**d*(τ) or *d*^{β/ν}*m*¯*d*(τ) against *d*^{1}^{/ν}*τ, 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*U** ^{d}*-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 > T*

*), but this is to be expected: we are seeing scaling violations (a satisfactory fit to a quartic can be obtained (χ*

^{c}^{2}

_{29−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 precision of our simulation analysis.

Sample realisations from the model, taken at temperatures around the crit- ical 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 Γ

*of polygon graphs with no boundary vertices. Let Θ*

_{D}*be the subspace of Γ*

_{D}

^{W}*of graphs made up of just one polygon. Each point in Θ*

_{D}*corresponds to a single polygon, lying wholly in*

_{D}*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

*A(B** ^{χ}*)

*≤*X

*θ**⊂**γ*(*χ*)

*L(θ)*^{2}

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* *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
*µ*^{B}_{D}^{|}^{W}*≤*

Z

Θ*D*

*L(θ)*^{2}

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

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

R

Γ^{W}*D**δ(θ⊂γ)×e*^{−2}^{L}^{(}^{γ}^{)}^{/T}*dλ(γ)*
R

Γ^{W}*D**e*^{−2}^{L}^{(}^{γ}^{)}^{/T}*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}*:

*γ⊃θ}*

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

˜Γ^{θ}*W**D**≡* [

*γ**∈Γ*^{θ}*W**D*

*{γ\θ}*

is the sub-domain of interest. We have Z

Γ^{W}*D*

*e*^{−2}^{L}^{(}^{γ}^{)}^{/T}*dλ(γ)≥*
Z

˜Γ^{θ}*W**D*

*e*^{−2}^{L}^{(}^{γ}^{)}^{/T}*dλ(γ)* (8)
We now turn to the numerator of Equation (7). Carrying out the integration
over vertices in*θ* using the*δ-function,*

Z

Γ^{W}*D*

*δ(θ⊂γ)×e*^{−2}^{L}^{(}^{γ}^{)}^{/T}*dλ(γ) =*
Z

Γ^{θ}*W**D*

*e*^{−2}^{L}^{(}^{γ}^{)}^{/T}*κ(θ)dλ(γ\θ)*

= *κ(θ)e*^{−2}^{L}^{(}^{θ}^{)}* ^{/T}*
Z

˜Γ^{θ}*W**D*

*e*^{−2}^{L}^{(}^{γ}^{)}^{/T}*dλ(γ),* (9)
since*γ*does not contain*θ* in the second line. Substituting with (8) and (9) in
(7), and cancelling,

E{*δ(θ⊂γ)|∂D ∩B** ^{χ}*=

*∅ } ≤κ(θ)e*

^{−2}

^{L}^{(}

^{θ}^{)}

^{/T}*,*and consequently,

*µ*^{B}*D*^{|}^{W}*≤* 1
4πA(D)

Z

Θ*D*

*L(θ)*^{2}*e*^{−2}^{L}^{(}^{θ}^{)}^{/T}*dλ(θ).*

In close analogy with Griffiths’ proof, we obtain
*µ*^{B}*D*^{|}^{W}*≤* 1

4πA(D)
Z _{∞}

0 *b*^{2}*e*^{−2}* ^{b/T}*
Z

Θ*D*

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

*db* (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*Lemma. Let Θ

^{(}

_{D}

^{n}^{)}be the subset of Θ

*of polygons with*

_{D}*n*vertices.

Lemma *Let*

*J*^{n}*≡*
Z

Θ^{(}_{D}^{n}^{)}*δ(b−L(θ))dλ(θ)* (11)

*so that* Z

Θ*D*

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

*n*=3

*J**n*

*in (10). Then*

*J**n**≤A(D)n*^{2}(n*−*1)(2π)^{n}^{−1}*b*^{n}* ^{−3}*
(n

*−*2)!

*,*

*and consequently*

Z

Θ*D*

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

Proof of the Lemma: Start with *J**n* defined in Equation (11). Use a
standard labelling with*x*_{1} 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*
at*x*_{1} to make a chain. Consider the set ˜Θ^{(}_{D}^{n}^{)}of distinct non-intersecting chains
*θ*˜of*n* edges linking *n*+ 1 vertices, labeled with variables *x*_{1} to *x*^{n}_{+1}. All the
vertices in a chain lie entirely to the right of the first vertex (or directly above).

Polygons are chains, Θ^{(}_{D}^{n}^{)} *⊂* Θ˜^{(}_{D}^{n}^{)}, since the first and last vertices in a chain
may coincide. Transform variables from *{x*^{i}*}*^{n}*i*=1 to *{x*1*,{e**i**}*^{n}*i*=1*}, where* *e**i* is
a Cartesian vector with origin *x** ^{i}* corresponding to the edge from the

*i’th to*the (i+ 1)’th vertex. When we switch to integrating over chains, we constrain

*e*

_{1}+

*e*

_{2}+

*. . .*+

*e*

*n*to be zero, so that the polygon closes. Equation (11) becomes

*J**n**≤*
Z

Θ˜^{(}_{D}^{n}^{)}*δ(b−L(˜θ))δ*^{(2)}(Σ*k**e**k*) *de*_{1}*de*_{2}*. . . de**n*

*e*_{1}*e*_{2}*. . . e**n* *dx*_{1}*,*
with*e**i**≡L(e**i*) and using sin(ψ*i*)*≤*1.

The integrand is unbounded. We partition the space into regions, and impose
the constraints*b*=*L(˜θ) ande*_{1}+*e*_{2}+*. . .*+*e**n*= 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 *de**i*, 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}^{n}*,*^{)}*−* = *{θ*˜*∈*Θ˜^{(}_{D}^{n}^{)} *|* abs(L(˜*θ)−b)≥}*

Θ˜^{(}_{D}^{n,i}*,*^{)} = *{θ*˜*∈*Θ˜^{(}_{D}^{n}^{)} *|* abs(L(˜*θ)−b)< , e*^{i}*≥*(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* *n*edges and a total length differing from*b*by
not more than *, must have such an edge. The sets ˜*Θ^{(}_{D}^{n,i}*,*^{)}, *i* = 1,2*. . . n* are
not disjoint, but combine with ˜Θ^{(}_{D}^{n}*,*^{)}*−* to cover ˜Θ^{(}_{D}^{n}^{)}. Chains in Θ^{(}_{D}^{n}*,*^{)}*−* will not

contribute to the integral. It follows that
*J**n* *≤*

X*n*

*i*=1

Z

Θ˜^{(}_{D}^{n,i}*,*^{)}

*δ(b−L(˜θ))δ*^{(2)}(Σ*k**e**k*) *de*_{1}*de*_{2}*. . . de**n*

*e*_{1}*e*_{2}*. . . e**n* *dx*_{1}

*≤* *n*

(b*−)*
X*n*

*i*=1

Z

Θ^{(}_{D}^{n,i}*,*^{)}

*δ(b−L(θ))de*_{1}*de*_{2}*. . . de*_{−}*i**. . . de**n*

*e*_{1}*e*_{2}*. . . e*_{−}^{i}*. . . e*^{n}*dx*_{1}*.* (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 *x**n*+1 = *x*_{1}). We have carried out the integral *de**i**δ*^{(2)}(Σ*k**e**k*) and
used the bound on*e**i*.

The third step is to eliminate an edge length parameter, using*b*=*L(˜θ), the*
length constraint. Let*φ** ^{i}* denote the angle made by edge

*e*

*i*to a fixed direction in the plane. In polar coordinates Equation (13) is

*J*^{n}*≤* *n*
(b*−)*

X*n*

*i*=1

Z

Θ^{(}_{D}^{n,i}*,*^{)}

*δ(b−L(θ))de*_{1}*dφ*_{1}*. . . de*_{−}^{i}*dφ*_{−}^{i}*. . . de*^{n}*dφ*^{n}*dx*_{1}*.* (14)
For the polygon to close

*e** ^{i}*sin(φ

*) =*

^{i}*−*X

*n*

*k*=1
*k**6=**i*

*e** ^{k}*sin(φ

*), (15)*

^{k}*e** ^{i}*cos(φ

*) =*

^{i}*−*X

*n*

*k*=1
*k**6=**i*

*e** ^{k}*cos(φ

*), (16)*

^{k}and consequently

*b−L(θ) =b−*
X*n*

*k*=1
*k**6=**i*

*e**k**−e**k*cos(φ*k**−φ**i*).

Integrating *de**j* 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 the*jth 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=i*are not disjoint but their
union covers Θ^{(}_{D}^{n,i}*,*^{)}. From Equation (14)

*J*^{n}*≤* *n*
(b*−)*

X*n*

*i*=1

X*n*

*j*=1
*j**6=**i*

Z

Θ^{(}_{D}^{n,i,j}*,* ^{)}

*δ(b−L(θ))de*_{1}*dφ*_{1}*. . . de*_{−}^{i}*dφ*_{−}^{i}*. . . de*^{n}*dφ*^{n}*dx*_{1}*.*

We may now apply the integral *de** ^{j}* to the delta-function

*δ(b−L(θ)). We*transform from

*e, φ*to

*e*

^{0}*, φ*

*where*

^{0}*φ*

^{0}*k*=

*φ*

*and*

^{k}*e*

^{0}*k*=

*e*

*for 1*

^{k}*≤k*

*≤n,*

*k6=j,*and

*φ*

^{0}*j*=

*φ*

*and*

^{j}*e*^{0}*j*=*e*^{j}*−e** ^{j}*cos(φ

^{j}*−φ*

*).*

^{i}The Jacobian of the full transformation*e, φ→e*^{0}*, φ** ^{0}* is just

*J*

*(e, φ*

^{−1}*→e*

^{0}*, φ*

*) =*

^{0}*∂e*

^{0}*j*

*∂e**j*

= 1*−*cos(φ^{j}*−φ** ^{i}*)

*−e*

*sin(φ*

^{j}

^{j}*−φ*

*)*

^{i}*∂φ*

^{i}*∂e**j**.* (17)
Repeated use of Equations (15) and (16) gives

*∂φ*^{i}

*∂e** ^{j}* =

*−*sin(φ

^{j}*−φ*

*)*

^{i}*e*^{i}*,*

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(θ))de** ^{j}*by one, and restrict the integration domain to polygons of length

*b,ie*set = 0. We obtain the simplified bound

*J*^{n}*≤* *n*
*b*

X*n*

*i*=1

X*n*

*j*=1
*j**6=**i*

Z

Θ^{(}_{D}^{n,i,j}*,*=0^{)}

*de*_{1}*dφ*_{1}*. . . de*_{−}^{j}*dφ*^{j}*. . . de*_{−}^{i}*dφ*_{−}^{i}*. . . de*^{n}*dφ*^{n}*dx*_{1}*.* (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*_{1}to remain in*D. This will*
include polygons with crossing edges and allow the polygon to overlap the border
of*D. The integraldx*_{1} 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

*e*_{1}+*e*_{2}+*. . . e*_{−}*j*+*. . . e*_{−}*i*+*. . .*+*e**n**≤b−b/n*

of volume less than *b*^{n}^{−2}*/(n−*2)!. Combining these factors with a factor of
(n*−*1) from the sum over*j, we obtain the bound on* *J**n* 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).

### 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.

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. *Proba-*
*bility Theory Related Fields, 80:543–579, 1989.*

[13] RB Griffiths. Peierls’ proof of spontaneous magnetization of a two dimen-
sional 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.*

[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).

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 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.

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

dd=6
d=8
d=12
d=16
T^{*},U^{*}

Figure 4: Binder parameter *U**d* (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σ.

-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

dFigure 5: The Binder parameter data of Figure 4 rescaled with Ising critical
exponents. The regression is a cubic polynomial. *χ*^{2}_{43−4} = 38.5 for the fit is
acceptable.

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

dFigure 6: The magnetisation ¯*m** ^{d}*(T), regressed with cubic polynomials.

-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

dFigure 7: The magnetisation data of Figure 6 rescaled with Ising critical ex-
ponents. The regression is a quartic polynomial. The value of the*χ*^{2} statistic
shows that the fit is a poor one.

### T=0.62 0.64

### 0.70

### 0.66

### 0.72 0.68

Figure 8: A selection of states equilibrated in a box of side*d*= 12 at tempera-
tures below and above the estimated critical temperature*T*^{c}*'*0.6665(5).