PERFECT SIMULATION OF THE HARD DISKS MODEL
BY PARTIAL REJECTION SAMPLING
HENG GUO AND MARK JERRUM

arXiv:1801.07342v1 [math.PR] 22 Jan 2018

Abstract. We present a perfect simulation of the hard disks model via the partial rejection
sampling method. Provided the density of disks is not too high, the method produces exact
samples in O(log n) rounds, where n is the expected number of disks. The method extends
easily to the hard spheres model in d > 2 dimensions.

1. Introduction
The hard disks model is one of the simplest gas models in statistical physics. Its configurations
are non-overlapping disks of uniform radius r in a bounded region of R2 . For convenience, in
this paper, we take this region to be the unit square [0, 1]2 . This model was precisely the one
studied by Metropolis et al. [11], in their pioneering work on the Markov chain Monte Carlo
(MCMC) method. They used Los Alamos’ MANIAC computer to simulate a system with 224
disks.
There are two variants of this model. To obtain the canonical ensemble, we fix the number
(or equivalently, density) of disks and decree that all configurations are “equally likely”, subject
only to the disks not overlapping. In the grand canonical ensemble, we fix the “average” number
of disks. To be more specific, centers of the disks are distributed according to a Poisson point
process of intensity λr = λ/(πr 2 ), conditioned on the disks being non-overlapping. The hard
disks model, and its higher dimensional generalization (called the hard spheres model) are also
related to the optimal sphere packing density [4, 17, 2]. See [6, 1] and references therein for
more details.
Our main aim in this work is to describe and analyse a very simple algorithm for exactly
sampling from the grand canonical ensemble, based on the partial rejection sampling paradigm
introduced by Guo, Jerrum and Liu [3].
More precisely, the challenge is the following: produce a realisation P ⊂ [0, 1]2 of a Poisson
point process of intensity λr in the unit square, conditioned on the event that no pair of points
in P are closer than 2r in Euclidean distance. We refer to this target measure as the hard
disks distribution. It describes an arrangement of open disks of radius r with centres in [0, 1]2
that are not allowed to overlap, but which otherwise do not interact. It is a special case of
the Strauss process [16]. Note that, although the disks do not overlap each other, they may
extend beyond the boundary of the unit square. Also, the intensity of the underlying Poisson
process is normalised so that the expected number of points of P lying in a disk of radius r
is 1. This normalisation gives us sensible asymptotics as the radius of the disks tends to zero
(equivalently, the number of disks tends to infinity).
Classical rejection sampling applied to this problem yields the following algorithm: repeatedly
sample a realisation P of the Poisson process of intensity λ in the unit square until P satisfies
the condition that no two points are closer than 2r, and return P . Unfortunately, for every
λ > 0, however small, the expected number of unsuccessful trials using this approach increases
exponentially in r −1 , as r → 0. Partial rejection sampling [3] requires only a subset of P to
be resampled at each iteration. Algorithm 1 below arises from a routine application of the
paradigm to the problem at hand.
The original partial rejection method [3] and its analysis are tailored for the discrete case. In
this paper we provide an alternative view on the correctness of the method, which is also valid
The work described here was supported by the EPSRC research grant EP/N004221/1 “Algorithms that
Count”.
1

2

HENG GUO AND MARK JERRUM

in the continuous setting. In other words, as with classical rejection sampling, Algorithm 1
terminates with probability 1, producing a realisation of the exact hard disks distribution. In
contrast to classical rejection sampling, the expected number of iterations (resampling steps)
is now asymptotically O(log(r −1 )) as r → 0, provided λ is not too large. We prove that rapid
termination occurs when λ < 0.21027. This analysis is not tight, and experiments suggests that
the actual threshold for rapid termination is around λ ≈ 0.51. The experimental advantage of
partial rejections is its simple termination rule, unlike MCMC, where it is difficult to determine
how long the algorithm should run.
The method extends naturally to the hard spheres model in d > 2 dimensions. Here, the
desired distribution is a Poisson point process in [0, 1]d conditioned on no pair of points being
closer than 2r. The natural normalisation for the intensity of the Poisson process in d dimensions
is λr,d = λ/(vd r d ), where vd is the volume of a ball of radius 1 in Rd . With this convention, we
1
prove that rapid termination occurs in d dimensions provided λ < 2−(d+ 2 ) .
The expected packing density α for this model is the expected total volume of balls. (Note that,
neglecting boundary effects, α is the proportion of the unit cube occupied by balls.) Intuitively,
we expect α to grow monotonically with λ, but that its rate of growth will slow dramatically as
the balls pack more tightly. The connection between expected packing density α and intensity λ
has recently been thoroughly explored by Jenssen, Joos and Perkins [6]. Using their results,
we show that partial rejection sampling can achieve expected packing density Ω(2−d ) while
retaining runtime O(log(r −1 )). Although sphere packings of density Ω(d2−d ) have been proved
to exist, there is no polynomial-time sampler that provably achieves packing density beyond
O(2−d ), as far as we are aware.
Other approaches to exact sampling include Coupling From The Past (CFTP), which was
adapted to point processes by Kendall [8] and Kendall and Møller [9]. Recently, Moka, Juneja
and Mandjes [12] proposed an algorithm based on rejection and importance sampling. Although
this algorithm, like ours, is based on rejection sampling, it does not share our asymptotic
performance guarantee. Indeed, its running time appears to grow exponentially as the number
of disks goes to infinity, with the density of disks held constant. See also [6, 13] for rigorous
bounds on packing densities of grand canonical ensembles.
Approximate sampling via Markov chain simulation has been studied by Kannan, Mahoney
and Montenegro [7] and Hayes and Moore [5] in the context of the canonical ensemble, where
the number of disks in a configuration is fixed. The best rigorous density bound to guarantee
rapid mixing is ≈ 0.154 [5]. It is should be noted that this is not directly comparable with our λ
due to the difference in models.
2. The sampling algorithm
The following notation will be used throughout. If P is a finite subset of [0, 1]2 then

BadPairs(P ) = {x, y} : x, y ∈ P ∧ x 6= y ∧ kx − yk < 2r ,

where k · k denotes Euclidean norm, and

BadPoints(P ) =

[

BadPairs(P ).

The open disk of radius r with centre x ∈ [0, 1]2 is denoted by Dr (x). The finite set Π ⊂ [0, 1]2
always denotes a realisation of the Poisson point process of intensity λr on [0, 1]2 . For a random
variable X and event E we use D(X) to denote the distribution (law) of X, and D(X | E) the
distribution of X conditioned on E occurring. Thus, D(Π | BadPoints(Π) = ∅) is the hard disks
distribution that we are interested in.
Our goal is to analyse the correctness and running time of a sampling algorithm for the hard
disks model (see Algorithm 1 below), specifically to determine the largest value of λ for which
it terminates quickly, i.e., in O(log r −1 ) iterations. The algorithm is an example application of
“Partial Rejection Sampling” [3], adapted to the continuous state space setting.

PERFECT SIMULATION OF THE HARD DISKS MODEL BY PARTIAL REJECTION SAMPLING

3

Algorithm 1 Partial Rejection Sampling for the hard disks model
PRS for Hard Disks(λ, r) // r is the disk radius and λr = λ/(πr 2 ) the intensity
Let P be a sample from the Poisson point process of intensity λr on the unit square
while B ← BadPoints(P ) 6= ∅ do
S ← B + D2r (0) // Resampling set is the Minkowski sum of B with a disk of radius 2r
Let PS be a sample from the Poisson point process of intensity λr on S
P ← (P \ B) ∪ PS
end while
return P
3. Correctness
Let B be any finite subset of [0, 1]2 . We say that B is a feasible set of bad points if
BadPoints(B) = B; this is equivalent to saying that there is a finite subset R ⊂ [0, 1]2 with
B = BadPoints(R). The key to establishing correctness of Algorithm 1 is the following loop
invariant:
D(P | BadPoints(P ) = B) = D(Π | BadPoints(Π) = B),
for every feasible set B, where P is any intermediate set of points during the execution of the
algorithm, and Π is a realisation of the Poisson point process. Let us consider what the right
hand side means operationally. Let S = B + D2r (0). (This is the resampling set used by the
algorithm.) Let Q be a sample from the distribution D(Π | BadPoints(Π) = B). The only points
in Q that lie inside S are the points in B. (Any extra points would create more bad pairs than
there actually are.) Thus Q ∩ S = B. Outside of S there are no bad pairs; thus Q ∩ S is a
sample from the hard disks distribution on S = [0, 1]2 \ S. Note that, setting B = ∅, we see
that D(Π | BadPoints(Π) = ∅) is just the hard disks distribution on [0, 1]2 .
Theorem 1. Algorithm 1 is correct: conditional on halting, Algorithm 1 produces a sample
from the hard disks distribution with intensity λr = λ/(πr 2 ).
Theorem 1 follows from Lemma 2 below.
Let T (a random variable) be the number of iterations of the while-loop. On each iteration, the
while loop terminates with probability bounded away from 0; thus T is finite with probability 1.
(Indeed, T has finite expectation.) Let Pt , for 1 ≤ t ≤ T , be the point set P after t iterations
of the loop, and P0 be the initial value of P (which is just a realisation of the Poisson point
process on [0, 1]2 ). We say that B0 , B1 , . . . , Bt ⊂ [0, 1]2 is a feasible sequence of (finite) point
sets if there exists a run of Algorithm 1 with BadPoints(P0 ) = B0 , . . . , BadPoints(Pt ) = Bt .
Lemma 2. Let B0 , B1 , . . . , Bt ⊂ [0, 1]2 be a feasible sequence. Then

D Pt BadPoints(P0 ) = B0 ∧ · · · ∧ BadPoints(Pt ) = Bt = D(Pt BadPoints(Pt ) = Bt )
= D(Π BadPoints(Π) = Bt ).

Proof. We prove the result by induction on t. The base case, t = 0, holds by construction: P0
is just a realisation of the Poisson point process on [0, 1]2 . Our induction hypothesis is

(1)
D Pt BadPoints(P0 ) = B0 ∧ · · · ∧ BadPoints(Pt ) = Bt = D(Π BadPoints(Π) = Bt ),

for every feasible sequence B0 , . . . , Bt . Extend the feasible sequence to Bt+1 . For the inductive
step, we assume (1) and aim to derive

(2) D Pt+1 BadPoints(P0 ) = B0 ∧ · · · ∧ BadPoints(Pt+1 ) = Bt+1
= D(Π | BadPoints(Π) = Bt+1 ).

The resampling set on iteration t + 1 is S = Bt + D2r (0). As a first step we argue below that

(2)
(3) D Pt+1 BadPoints(P0 ) = B0 ∧ · · · ∧ BadPoints(Pt ) = Bt = D(Π | BadPairs(Π) ∩ S = ∅),
(2)

where S
denotes the set of unordered pairs of elements from S. We have noted that (1)
implies that, outside of the resampling set S, the point set Pt is a realisation of the hard disks

4

HENG GUO AND MARK JERRUM

distribution. Also, the algorithm does not resample points outside of S. Thus Pt+1 ∩ S = Pt ∩ S
is Poisson distributed, conditioned on there being no bad pairs. Inside S, resampling has left
behind a fresh Poisson point process Pt+1 ∩ S. These considerations give (3).
Next, we condition on Bt+1 :

D Pt+1 BadPoints(P0 ) = B0 ∧ · · · ∧ BadPoints(Pt ) = Bt ∧ BadPoints(Pt+1 ) = Bt+1

(2)
= D Π BadPairs(Π) ∩ S = ∅ ∧ BadPoints(Π) = Bt+1 .

Since Bt+1 contains no bad pairs with both endpoints in S, the event BadPoints(Π) = Bt+1
(2)
entails the event BadPairs(Π) ∩ S = ∅. Thus, we have

D Pt+1 BadPoints(P0 ) = B0 ∧ · · · ∧ BadPoints(Pt ) = Bt ∧ BadPoints(Pt+1 ) = Bt+1
= D(Π | BadPoints(Pt+1 ) = Bt+1 ).

The right hand side of this equation does not involve B0 , . . . , Bt , and so
D(Pt+1 | BadPoints(Pt+1 ) = Bt+1 ) = D(Π | BadPoints(Π) = Bt+1 ).
This completes the induction step (2) and the proof.



Proof of Theorem 1. As we observed earlier, T , the number of iterations of the while-loop, is
finite with probability 1. By Lemma 2, noting that BT = ∅,
D(PT ) = D(Π | BadPoints(Π) = ∅).

In other words, at termination, Algorithm 1 produces a realisation of the hard disks distribution
on [0, 1]2 .

4. Run-time analysis
We consider how the number of “bad events”, i.e., the cardinality of the set BadPairs(Pt ),
evolves with time. As usual Π denotes a realisation of the Poisson point process of intensity λr .
Also denote by ∆ a realisation of the hard disks process of the same intensity. We need the
following stochastic domination result.
Lemma 3. The hard disks distribution is stochastically dominated by the Poisson point process
with the same intensity. That is, we can construct a joint sample space for Π and ∆ such that
∆ ⊆ Π.
Holley’s criterion is a useful test for stochastic domination, but it is not of direct use to us in
the proof of Lemma 3, because it applies only to finite state spaces. Fortunately, Preston [14,
Theorem 9.1], has derived a version of Holley’s criterion that fits our situation. We will mostly
follow Preston’s notation, except that, to save confusion, we will use P and Q, rather than x
and y, to denote finite sets of points. In order to state his result, we need some notations. In our
application, ω
en is the distribution on ([0, 1]2 )(n) obtained by sampling n points independently
and uniformly
at random from [0, 1]2 , and regarding the points as indistinguishable; furthermore,
P∞
ω
e = n=0 ω
en /n!. Denote by Ω the set of all finite subsets of [0, 1]2 , and by F the set of nonnegative measurable functions Ω → R satisfying
Z
(4)
f de
ω = 1,
and
(5)

f (P ) = 0 and P ⊆ Q implies f (Q) = 0.

(See Preston [14, Section 9] for detailed formal definitions of the concepts here.)
Lemma 4 (Theorem 9.1 of [14]). Let f1 , f2 ∈ F and suppose that
(6)

f2 (Q + ξ)
f1 (P + ξ)
≥
,
f1 (P )
f2 (Q)

for all Q ⊆ P ∈ Ω and ξ ∈ [0, 1]2 \ P

PERFECT SIMULATION OF THE HARD DISKS MODEL BY PARTIAL REJECTION SAMPLING

5

(where, by convention, 0/0 = 0). Then, for all bounded, measurable, non-decreasing functions
g : Ω → R,
Z
Z
gf1 de
ω ≥ gf2 de
ω,
i.e., if µi is the probability measure having density fi with respect to ω
e , then µ1 stochastically
dominates µ2 .
Proof of Lemma 3. We set

|
f1 (P ) = C1 λ|P
r

and

(
|P |
C 2 λr ,
f2 (P ) =
0,

if BadPairs(P ) = ∅;
otherwise.

The normalising constants C1 and C2 are chosen so that both f = f1 and f = f2 satisfy (4).
(There is an explicit expression for C1 , namely C1 = exp(−λr ), but not for C2 .) Notice that
both f1 and f2 also satisfy (5). Notice also that the probability measures µ1 and µ2 of the
Poisson point process and the hard disks process have densities f1 and f2 with respect to ω
e.
The premise (6) of Lemma 4 holds, since the left-hand side is always λr and the right-hand side
is either λr or 0. The conclusion is that µ1 dominates µ2 . Strassen’s Theorem [10, 15], allows
us to conclude the existence of a coupling of Π and ∆ as advertised in the statement of the
lemma (except, possibly, on a set of measure zero).

Lemma 5. There is a bound λ̄ > 0 such that the expected number of iterations of the while-loop
in Algorithm 1 is O(log r −1 ) when λ < λ̄.
Proof. First observe that BadPairs(P ) determines BadPoints(P ) and vice versa. So conditioning
on the set BadPairs(P ) is equivalent to conditioning on BadPoints(P ).
Introduce random variables Zt = |BadPairs(Pt )|, for 1 ≤ t ≤ T . Our strategy is to show that
(7)

E(Zt+1 | Z0 , . . . , Zt ) ≤ α−1 Zt ,

for some α > 1. Then Z0 , αZ1 , α2 Z2 , α3 Z3 , . . . is a supermartingale (with the convention that
Zt = 0 for all t > T ). Therefore, E Zt ≤ α−t E Z0 ≤ 21 λ2r α−t . Here, we have used the fact that
|P0 | is a Poisson random variable with expectation λr , and
Z0 = |BadPairs(P0 )| ≤ 12 |P0 | (|P0 | − 1),

and hence



E Z0 ≤ E |P0 | (|P0 | − 1) = 12 λ2r .

Setting t = O(log r −1 + log ε−1 ), we obtain E Zt < 1/ε, and hence Pr(T > t) ≤ ε. It follows that
the expected number of iterations of the while-loop of Algorithm 1 is O(log r −1 ). Note that
the probability of non-termination decreases exponentially with t, so the probability of large
deviations above the expected value
√ of T is low.
Crude estimates give λ̄ = 1/(4 2). The calculation goes as follows. Suppose, in (7), we
condition on the random variables BadPoints(P0 ), . . . , BadPoints(Pt ), rather than simply on
Z0 , . . . , Zt . This is more stringent conditioning, since the former random variables determine the latter. It is enough to establish (7) under the more stringent conditioning. So fix
BadPoints(P0 ) = B0 , . . . , BadPoints(Pt ) = Bt , and note that this choice also fixes the resampling sets S0 , . . . , St . Suppose Zt = |BadPairs(Pt )| = k. Inside the resampling set St we have
a Poisson point process Pt+1 ∩ St of intensity λr . Outside, by Lemma 2, there is a realisation
of the hard disks process. Since we are seeking an upper bound on Zt+1 we may, by Lemma 3,
replace Pt+1 ∩ S t by a Poisson point process of intensity λr .
Let k′ = E(Zt+1 | Z0 , . . . , Zt ). From the above considerations we have
Z
Z
′
λr 1kx−yk≤2r dy dx.
λr
(8)
k ≤
St

[0,1]2

6

HENG GUO AND MARK JERRUM

This is an overestimate, as we are double-counting overlapping disks whose centres both lie
within St . Now, St is a union of at most 2k disks of radius 2r. Thus
(9)

′

k ≤

2kλ2r

=

2kλ2r

Z

Z

D2r (0)

Z

R2

Z

1kx−yk≤2r dy dx
dy dx

D2r (0) D2r (x)
2
2

= 2kλ2r × 4πr × 4πr
= 32λ2 k.

There are further sources of slack here: there may be fewer than 2k disks, the disks comprising
St certainly overlap, and, for points x near the boundary, some of disks D2r (x) will lie partly
outside the unit square.
(The last of these presumably has no effect asymptotically, as r → 0.)
√
Setting λ̄ = 1/(4 2) = 0.17677+, we see that α = k/k′ > 1 for any λ < λ̄, and (Zt )∞
t=0 is a
supermartingale, as required.

The constant λ̄ may seem quite small. Note, however, that classical rejection sampling cannot
achieve any λ̄ > 0. The argument goes as follows. Divide [0, 1]2 into r × r squares. If there
are two points in the same square then they will certainly be less than distance 2r apart. The
number of points in each square is Poisson distributed with mean λ/π. Thus for any λ > 0 the
probability that a particular square has at least two points is bounded away from zero. The
number of points in each square is independent of all the others. It follows that the running
time of classical rejection sampling is exponential in r −2 .
The above derivation for λ̄ is quite crude and can be improved.
Lemma 6. The constant λ̄ in Lemma 5 can be taken to be 0.21027.
Proof. For each of the 2k disks, the right-hand side of inequality (9) counts pairs of points (x, y)
with x in the disk, and y anywhere within distance 2r of x. Since a bad event is determined by
an unordered pair of points, this gives rise to significant double counting. In particular, pairs
(x, y) with x and y lying in the same ball are double counted. We can subtract off these pairs
to obtain a better estimate.
For a single ball, the correction is
1
C=
2

Z

λr

D2r (0)

Z

D2r (0)

λr 1kx−yk≤2r dy dx.

(The initial factor of one half arises because we want to count unordered pairs.) With the
change of variables x = 2rx′ and y = 2ry ′ this expression simplifies to
Z
Z
1
4 2
1kx′ −y′ k≤1 dy ′ dx′
C = × 16r λr
2
D1 (0) D1 (0)
Z
= 8λ2r r 4
L(kx′ k) dx′ ,
D1 (0)

where L(kx′ k) is the area of the “lens” D1 (0) ∩ D1 (x′ ). Letting ̺ denote the offset of the centres
of the two disks, the area of the lens is given by
L(̺) = 2 arccos(̺/2) − 12 ̺

p

4 − ̺2 .

(This is by elementary geometry: the lens is the intersection of two sectors, one from each of the
disks, and its area can be computed by inclusion-exclusion.) Translating to polar coordinates

PERFECT SIMULATION OF THE HARD DISKS MODEL BY PARTIAL REJECTION SAMPLING

(̺, θ),
C=

8λ2r r 4

= 8λ2r r 4

Z

7

L(kx′ k2 ) dx′

D1 (0)
Z 2πZ 1
0

̺ L(̺) d̺ dθ

0
Z 1

8λ2
̺ L(̺) d̺
× 2π
π2
0
1

 ̺   ̺ ̺3 p
16λ2 π
2
2
4−̺
=
−
+ (̺ − 1) arccos
+
π
2
2
4
8
0
√ 

6 3
= 8−
λ2 .
π

=

(The integral was evaluated using the Maple computer algebra system.) Our revised upper
bound on k′ is thus
√ 

6 3
′
2
2
k ≤ 2k(16λ − C) = 2kλ 8 +
(10)
,
π
Solving

√ 
12 3
=1
λ̄ 16 +
π
yields the improved bound of λ̄ = 0.21027+.
2





There are other factors that could in principle be used to increase λ̄ further — each disk
necessarily overlaps with at least one other disk, some bad events are triple or quadruple counted
— but the computational difficulties rapidly increase when attempting to account for these.
5. Three or more dimensions
In higher dimensions, the hard disk model is known as the hard spheres model. Everything
in Sections 3 and 4 carries across to d > 2 dimensions with little change. For general d, the
appropriate scaling for the intensity is λr,d = λ/(vd r d ), where vd is the volume of a ball of unit
radius in d dimensions. Note that in a realisation of a Poisson point process with intensity λr,d ,
the expected number of points in a ball of radius r is λ.
The analogue of equation (8) is
Z
Z
′
λr,d 1kx−yk≤2r dy dx,
λr,d
k ≤
St

[0,1]d

which leads to

k′ ≤ 22d+1 λ2 k.

1

So setting λ̄ = 2−(d+ 2 ) we find that α = k/k′ > 1 for any λ < λ̄. It follows that the running
time of partial rejection sampling is O(log r) for any λ < λ̄.
By a result of Jenssen, Joos and Perkins [6], we lose just a constant factor when translating
from intensity λ to packing density α. (It is partly to connect with their work, we measure
intensity in terms of the expected number of points in a ball of radius r.) In the proof of [6,
Thm 2], the following inequality is derived:

α ≥ inf max λe−z , 2−d exp[−2 · 3d/2 λ] · z .
z
√
Assuming λ ≤ λ̄, which holds in the range of validity of our algorithm, we have 2λ ≤ 2−d and
hence
√
√

α ≥ inf max λe−z , 2λ exp[− 2(3/4)d/2 ] · z
z

= cd λ,

where

8

HENG GUO AND MARK JERRUM

√
√

cd = inf max e−z , 2 exp[− 2(3/4)d/2 ] · z .
z

Note that (cd ) is monotonically increasing, with c2 = 0.42220+, and limd→∞ cd = 0.63724+. It
follows that we can reach expected packing density Ω(2−d ) with O(log r −1 ) expected iterations.
This is currently the best that can be achieved by any provably correct sampling algorithm with
polynomial (in 1/r) runtime [7]. The asymptotically best packing density currently rigorously
known is d2−d , but achieving this would require λ to grow exponentially fast in d. This is
clearly beyond the capability of partial rejection sampling, but also beyond the capability of
any known efficient sampling algorithm.
Acknowledgements
We thank Mark Huber and Will Perkins for inspiring conversations and bringing the hard
disks model to our attention.
References
[1] Henry Cohn. A conceptual breakthrough in sphere packing. Notices Amer. Math. Soc., 64(2):102–115, 2017.
1
[2] Henry Cohn, Abhinav Kumar, Stephen D. Miller, Danylo Radchenko, and Maryna Viazovska. The sphere
packing problem in dimension 24. Ann. of Math. (2), 185(3):1017–1033, 2017. 1
[3] Heng Guo, Mark Jerrum, and Jingcheng Liu. Uniform sampling through the Lovasz local lemma. In STOC,
pages 342–355, 2017. 1, 2
[4] Thomas C. Hales. A proof of the Kepler conjecture. Ann. of Math. (2), 162(3):1065–1185, 2005. 1
[5] Thomas P. Hayes and Cristopher Moore. Lower bounds on the critical density in the hard disk model via
optimized metrics. CoRR, abs/1407.1930, 2014. 2
[6] Matthew Jenssen, Felix Joos, and Will Perkins. On the hard sphere model and sphere packings in high
dimensions. ArXiv, abs/1707.00476, 2017. 1, 2, 7
[7] Ravi Kannan, Michael W. Mahoney, and Ravi Montenegro. Rapid mixing of several Markov chains for a
hard-core model. In ISAAC, pages 663–675, 2003. 2, 8
[8] Wilfrid S. Kendall. Perfect simulation for the area-interaction point process. In Probability towards 2000
(New York, 1995), volume 128 of Lect. Notes Stat., pages 218–234. Springer, New York, 1998. 2
[9] Wilfrid S. Kendall and Jesper Møller. Perfect simulation using dominating processes on ordered spaces, with
application to locally stable point processes. Adv. Appl. Probab., 32(3):844–865, 2000. 2
[10] Torgny Lindvall. On Strassen’s theorem on stochastic domination. Electron. Comm. Probab., 4:51–59, 1999.
5
[11] Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward
Teller. Equation of state calculations by fast computing machines. J. Chem. Phys., 21(6):1087–1092, 1953. 1
[12] S. B. Moka, S. Juneja, and M. R. H. Mandjes. Perfect sampling for Gibbs processes with a focus on hardsphere models. ArXiv, abs/1705.00142, 2017. 2
[13] Will Perkins. Birthday inequalities, repulsion, and hard spheres. Proc. Amer. Math. Soc., 144(6):2635–2649,
2016. 2
[14] Chris Preston. Spatial birth-and-death processes (with discussion). Bull. Inst. Internat. Statist., 46(2):371–
391, 405–408 (1975), 1975. 4
[15] V. Strassen. The existence of probability measures with given marginals. Ann. Math. Statist., 36:423–439,
1965. 5
[16] David J. Strauss. A model for clustering. Biometrika, 62(2):467–475, 1975. 1
[17] Maryna S. Viazovska. The sphere packing problem in dimension 8. Ann. of Math. (2), 185(3):991–1015,
2017. 1
(Heng Guo) School of Informatics, University of Edinburgh, Informatics Forum, Edinburgh, EH8
9AB, United Kingdom.
E-mail address: hguo@inf.ed.ac.uk
(Mark Jerrum) School of Mathematical Sciences, Queen Mary, University of London, Mile End
Road, London, E1 4NS, United Kingdom.
E-mail address: m.jerrum@qmul.ac.uk

