Sparse Hard-Disk Packings and Local Markov Chains

We propose locally stable sparse hard-disk packings, as introduced by Böröczky, as a model for the analysis and benchmarking of Markov-chain Monte Carlo (MCMC) algorithms. We first generate such Böröczky packings in a square box with periodic boundary conditions and analyze their properties. We then study how local MCMC algorithms, namely the Metropolis algorithm and several versions of event-chain Monte Carlo (ECMC), escape from configurations that are obtained from the packings by slightly reducing all disk radii by a relaxation parameter. We obtain two classes of ECMC, one in which the escape time varies algebraically with the relaxation parameter (as for the local Metropolis algorithm) and another in which the escape time scales as the logarithm of the relaxation parameter. A scaling analysis is confirmed by simulation results. We discuss the connectivity of the hard-disk sample space, the ergodicity of local MCMC algorithms, as well as the meaning of packings in the context of the NPT ensemble. Our work is accompanied by open-source, arbitrary-precision software for Böröczky packings (in Python) and for straight, reflective, forward, and Newtonian ECMC (in Go).

arbitrary-precision software for Böröczky packings and for ECMC. We discuss the apparent paradox that Böröczky packings, on the one hand, render local MCMC non-irreducible (that is, "non-ergodic") but on the other hand do not invalidate their practical use. We resolve this paradox by considering the N P T ensemble (where the pressure is conserved instead of the volume). We moreover advocate the usefulness of ε-relaxed Böröczky configurations for modeling bottlenecks in MCMC and propose the comparison of escape times from these configurations as a useful benchmark for real-world problems.
This work is organized as follows. In Section 2, we construct Böröczky packings following the original proposal [11] and a variant due to Kahle [15], and we analyze their properties. In Section 3, we discuss local MCMC algorithms and present analytical and numerical results for the escape times from the ε-relaxed Böröczky configurations. In Section 4, we analyze algorithms and their escape times and discuss fundamental aspects, among them irreducibility, statistical ensembles, the question of bottlenecks, and the difference between local and non-local MCMC methods. In the conclusion (Section 5), we point to several extensions and place our findings into the wider context of equilibrium statistical mechanics, the physics of glasses and the mechanics of granular materials. Our open-source arbitrary-precision software for Böröczky packings and for ECMC is presented in Appendix A.

Böröczky packings
In the present section, we consider packings of N disks of radius σ = 1 in a periodic square box of sides L. The density η is the ratio of the disk areas to that of the box: For concreteness, the central simulation box ranges from −L/2 to L/2 in both the x and the y direction. The periodic boundary conditions map the central simulation box onto an infinite hard-disk system with periodically repeated boxes or, equivalently, onto a torus. In a locally stable packing, each disk is blocked-at a distance 2σ-by at least three other disks (taking into account periodic boundary conditions), with the contacts not all in the same halfplane. The opening angle of a disk i, the largest angle formed by the contacts to its neighbors, is then always smaller than π. The maximum opening angle is the largest of the N opening angle of all disks. Clearly, the packing cannot be escaped from through the infinitesimal single-disk moves of ECMC or, in Metropolis MCMC, through steps of small enough range (see Subsection 3.1).

Construction of Böröczky packings
In the central simulation box, a finite-N Böröczky packing is built on a central core placed around (0, 0). This core connects to four of its periodic copies  Cores with different shapes, as for example that of a triangle, yield Böröczky packings in other geometries [11,16,15].

Branches-infinite-layer case (infinite N )
Following Ref. [11], we first construct infinite branches that correspond to the N → ∞ and η → 0 limits, without periodic boundary conditions. One such branch is attached to each of the four sides of the central core so that all disks are locally stable. The horizontal branch that extends from the central core in the positive x-direction is symmetric about the x-axis. The half branch for y ≥ 0 uses three sets of disks For the branch that is symmetric about the x-axis, the construction relies on four horizontal lines [11]: The disks A 1 and B 1 are aligned in x at heights g 3 and g 1 , respectively. All A disks lie on a given convex polygonal chain A between g 2 and g 3 . The chain segments on A are of length 2 so that subsequent disks A i and A i+1 block each other, and the position of A 1 fixes all other A disks. All C disks lie on g, and C i blocks B i from the right (in particular, C 1 is placed after B 1 ). The disk B i , for i > 1, lies between g and g 1 and it blocks disks A i and C i−1 from the right. With the position of g 2 , the branch approaches a hexagonal packing for i → ∞. After reflection about the x-axis, all disks except A 1 and B 1 are locally stable in the infinite branch. The Böröczky packing is completed by attaching the four branches along the four coordinate axes. For the Böröczky core, both A 1 and B 1 are blocked by core disks (see Fig. 1a). For the Kahle core, B 1 is blocked by a core disk, and A 1 is locally stable as it also belongs to another branch (see Fig. 1b).

Branches-finite-layer case (finite N ), periodic boundary conditions
Branches can also be constructed for periodic simulation boxes, with a finite number k of layers and finite N (see [11]). The branch that connects the central core placed around (0, 0) with its periodic image around (L, 0) is then again symmetric about the x-axis but, in addition, also about the boundary of the central simulation box at x = L/2. We describe the construction of the half-branch (for y ≥ 0) up to this boundary (see Fig. 1).
For half-branches with a finite number of layers k and a finite number of disks {A 1 , . . . , A k }, {B 1 , . . . , B k }, and {C 1 , . . . , C k−1 } (with their corresponding mirror images), the convex polygonal chain A lies between g < 2 and g 3 where g < 2 is an auxiliary horizontal line placed slightly below g 2 (see [11]). The horizontal lines g and g 1 and the algorithm for placing the disks are as in Subsection 2.1.2 (see Fig. 1c and d). By varying the distance between g 2 and g < 2 , one can make disk B k satisfy the additional requirement The position of B k then fixes the boundary of the square box (x B k = L/2) and B k blocks A k as well as the mirror image A k+1 of A k (see Fig. 1c

Properties of Böröczky packings
The BigBoro software package (see Appendix A) implements two convex polygonal chains that we now discuss. It also determines the collective escape modes from a Böröczky packing, the space of which we also discuss.

Convex polygonal chains (geometric, circular)
In the convex geometric chain A geo , the disks A i approach the line g < 2 exponentially in i. In contrast, in the circular chain A circ , all A disks lie on a circle (including their mirror images after reflection about x = L/2) so that their opening angles are all the same.
For the convex polygonal chain A geo , the distance between A i and g < 2 follows a geometric progression: with the attenuation parameter φ. (For a horizontal branch, the distances in eq. (3) are simply the difference between y-values.) The densities η Bör and η Kahle vary with φ, and they decrease as ∼ 1/k for large k (see Table 1). The geometric sequence for A i induces that the largest opening angle, usually the one between A k−1 , A k , A k+1 , approaches the angle π as θ k = φ k−2 (1 − φ)(g 3 − g < 2 )/2 ∼ φ k , that is, exponentially in k and in L.
For the convex polygonal chain A circ , all A disks lie on a circle of radius R, and in particular A 1 , which by construction is on g 3 (see Subsection 2.1.2). The circle is tangent to g < 2 at x = L/2. The center of the circle lies on the vertical line at x = L/2. It follows from basic trigonometry that for large k, the radius of the circle R scales as ∼ k 2 and the opening angles approach the angle π as ∼ k −2 .

Contact graphs: local and collective stability
The contact graph of a Böröczky packing connects any two disks whose pair distance equals 2 (possibly accounting for periodic boundary conditions, see Fig. 1). In a Böröczky packing with k ≥ 1 layers, the number N of disks and the number N contact of contacts are as follows: N N contact Böröczky core 20k + 12 32k + 20 Kahle core 20k − 4 32k + 4 .
For all values of k ≥ 1, the number of contacts is smaller than 2N − 2. This implies that collective infinitesimal two-dimensional displacements, with 2N − 2 degrees of freedom (the values of the displacements in x and in y for each disk avoiding trivial translations), can escape from a Böröczky packing, which is thus not collectively stable (see for example [15]). When all disks i, at positions x i , are moved to x i + ∆ i , the squared separation between disks i and j changes from |x i − x j | 2 to (possibly accounting for periodic boundary conditions). If the first-order term in eq. (5) vanishes for all i and j, the separation between disks in contact cannot decrease. It then increases-to second order in the displacements ∆ i , ∆ j -if ∆ i = ∆ j , so that contact is lost. The first-order term writes as a product of twice an "escape matrix" The row r of M esc corresponding to the contact between i and j has the following four non-zero entries The BigBoro software package (see Appendix A) solves for using singular-value decomposition. For the k = 5 Böröczky packing with the Kahle core, we find 28 vanishing singular values, that comprise the two uniform translations. It follows from eq. (4) that, because of 28 = 2N − N contact , all contacts are linearly independent. The corresponding space of all collective escape modes is 28-dimensional (see Fig. 2).
For all values of k ≥ 1, the number of contacts in eq. (4) is also larger than N − 1. Böröczky packings are thus collectively stable for displacements ∆ i that are constrained to a single direction (as for example the x or y direction). This impacts the behavior of unidirectional MCMC algorithms (see

Dimension of the space of Böröczky packings
As discussed in Subsection 2.2.2, each Böröczky packing has a contact graph. Conversely, a given contact graph may describe Böröczky packings for a continuous range of densities η. As an example, changing the attenuation parameter φ of the convex polygonal chain A geo in eq. (3) continuously moves all branch disks, and in particular disk B k and, therefore, the value of L and the density η (see Table 1 for density windows that can be obtained in this way). We expect that locally stable packings exist for any density at large enough N .
Moreover, the space B of locally stable packings of N disks of radius σ in a given central simulation box is of lower dimension than the sample space Ω: For each contact graph, each independent edge decreases the dimensionality by one. In addition there is only a finite number of contact graphs for a given N . The low dimension of B also checks with the fact that any packing, and more generally, any configuration with contacts, has effectively infinite pressure (see the detailed discussion in Subsection 4.2.2). As the ensemble-averaged pressure is finite (except for the densest packing), the packings must be of lower dimension. We conjecture Ω \ B to be connected for a given η below the densest packing at large enough N although, in our understanding, this is proven only for η ∼ 1/ √ N (see [17,18]).

MCMC algorithms and ε-relaxed Böröczky configurations
In this section, we first introduce to a number of local MCMC algorithms (see Subsection 3.1). In Subsection 3.2, we then determine the escape times (in the number of trials or events) after which these algorithms escape from ε-relaxed Böröczky configurations, that is, from Böröczky packings with disk radii multiplied by a factor (1 − ε). A scaling theory establishes the existence of two classes of MCMC algorithms, one in which the escape time from an ε-relaxed Böröczky configuration scales algebraically with ε, and the other in which the scaling is logarithmic. Numerical simulations confirm the theory.

Local hard-disk MCMC algorithms
We define the reversible Metropolis algorithm with two displacement sets, from which the trial moves are uniformly sampled. We also consider variants of the non-reversible ECMC algorithm that only differ in their treatment of events, that is, of disk collisions. An arbitrary-precision implementation of the discussd ECMC algorithms (in the Go programming language) is contained in the BigBoro software package (see Appendix A).

Local Metropolis algorithm: displacement sets
The N disks are at a position x = {x 1 , . . . , x N }. In the local Metropolis algorithm [1], at each time t = 1, 2, . . . , a trial move is proposed for a randomly chosen disk i, from its position x i to x i + ∆x i . If the trial produces an overlap, disk i stays put and x remains unchanged. We study two sets for the trial moves. For the cross-shaped displacement set, the trial moves are uniformly sampled within a range δ along the coordinate axes, that is, either along the x-axis (∆x i = (ran(−δ, δ) , 0)) or along the y-axis (∆x i = (0, ran(−δ, δ))).
Alternatively, for the square-shaped displacement set, the trial moves are uniformly sampled as ∆x i = (ran(−δ, δ) , ran(−δ, δ)). A Böröczky packing is invariant under the local Metropolis algorithm if the range δ is smaller than a critical range δ c . The latter is closely related to the critical opening angle (see the discussion in Subsection 2.2.1 and Fig. 3c). For these packings, the critical range vanishes for N → ∞. On the other hand, for large ranges δ, the algorithm can readily escape from the stable configuration. For δ = L/2, the Metropolis algorithm with a square-shaped displacement set proposes a random placement of the disk i inside the central simulation box. This displacement set leads to a very inefficient algorithm at the densities of physical interest, but it mixes very fast for sparse systems (see also Subsection 4.2.4).

Hard-disk ECMC: straight, reflective, forward, Newtonian
Straight ECMC [4] is one of the two original variants of event-chain Monte Carlo. This Markov chain evolves in (real-valued) continuous time t MCMC , but its implementation is event-driven. The algorithm is organized in a sequence of "chains", each with a chain time τ chain , its intrinsic parameter. During each chain time, disks move with unit velocity in one given direction (alternatively in +x or +y). A randomly sampled initial disk thus moves either until the chain time τ chain is used up, or until, at a collision event, it collides with another disk, which then moves in its turn, etc. This algorithm is highly efficient in some applications [4,5,19]. During each chain (in between changes of direction), any disk can collide only with three other disks or fewer [20,21]. A constraint graph with directed edges may encode these relations. This constraint graph (defined for hard-disk configurations) takes over the role of the contact graph (defined for packings) (see Fig. 3a and b). As the moves in a chain are all in the same direction, the straight ECMC has only N − 1 degrees of freedom, fewer than there are edges in the constraint graph. It is for this reason that it may encounter the rigidity problems evoked in Subsection 2.2.2. In reflective ECMC [4], in between events, disks move in straight lines just as in straight ECMC. At a collision event, the target disk does not continue in the same direction as the active disk. Rather, the target-disk direction is the original active-disk direction reflected from the line connecting the two disk centers at contact (see [4]). As all ECMC variants, reflective ECMC satisfies the global-balance condition. Because the number of disks is large, we need not implement resamplings as is necessary in low-dimensional systems in order to enforce irreducibility [22,23]. In consequence, the reflective ECMC has no intrinsic parameter. A variant of reflective ECMC, obtuse ECMC [14], has shown interesting behavior.
Forward ECMC [13], of which we implement the "Forward All Ref" variant, is a family of ECMC algorithms. After an event, the target-disk direction (of unit absolute value) is updated as follows. The component orthogonal to the line connecting the disks at contact is uniformly sampled between 0 and 1 (reflecting the orthogonal orientation). Its parallel component is determined so that the direction vector (which is also the velocity vector) is of unit norm. The parallel orientation remains unchanged. The forward ECMC has no intrinsic parameter and requires no resamplings.
Newtonian ECMC [14] mimics molecular dynamics in order to determine the velocity of the target disk in an event. It initially samples disk velocities from the two-dimensional Maxwell distribution. However, at each moment, only a single disk is actually moving with constant velocity. At a collision event, the velocities of the colliding disks are updated according to Newton's law of elastic collisions for hard disks of equal masses, but only the target disks actually moves after the event. In this algorithm, the velocity (which indexes the Monte-Carlo time) generally differs from unity. We do not implement resamplings, although (like reflective ECMC) Newtonian ECMC is not always irreducible without them [23]. As in earlier studies for three-dimensional hard-sphere systems [14], Newtonian ECMC is typically very fast for ε-relaxed Böröczky configurations. However, it suffers from frequent gridlocks (see Subsection 4.1.3).

Escape times from ε-relaxed Böröczky configurations
The principal figure of merit for a Markov chain is its mixing time [24], the number of steps it takes from the worst-case initial condition to approach the stationary probability distribution to some precision level. Böröczky packings are invariant under local Metropolis dynamics (of sufficiently small range) as well as under ECMC dynamics, so that the mixing times are, strictly speaking, infinite. Although they cannot be escaped from, the packings make up only a set of measure zero in sample space, and might thus be judged irrelevant.
However, as we will discuss in the present section, the situation is more complex. For every Böröczky packing, an associated ε-relaxed Böröczky configuration keeps the central simulation box and the disk positions, but reduces the disk radii from 1 to 1 − ε. An ε-relaxed Böröczky configuration effectively defines a finite portion of configuration space (the spheres of radius ε around each disk position, see Section 4.1.2). All MCMC algorithms considered in this work escape from these configurations in an escape times that diverges as ε → 0. We suggest that escape times are analogous to mixing times. In consequence, a finite portion of sample space is confined on times larger than an arbitrary constant. We suggest that the substantial differences between escape times may be relevant for real-world applications. The divergence of escape times for ε → 0 is specific to the N V T ensemble (see Subsection 4.2.2).

Nearest-neighbor distances and escape times
In a Böröczky packing, disks are locally stable, and they all have a nearestneighbor distance of 2. The packings are sparse, and the nearest-neighbor distance is thus smaller than its ∼ 1/ √ η equilibrium value. To track the escape from an ε-relaxed Böröczky configuration, we monitor the maximum nearestneighbor distance: where |x ij (t)| = |x j (t) − x i (t)| is the distance between disks i and j (possibly corrected for periodic boundary conditions). For the Metropolis algorithm, we compute d(t) once every N trials, and t denotes the integer-valued number of individual trial moves. For ECMC, we sample d(t) and the number of events in intervals of the sampling time. In eq. (8), t then denotes the integer-valued number of events. Starting from an ε-relaxed Böröczky configuration, d(t) typically remains at d(t) ∼ 2 + O (ε) for a long time until it approaches the equilibrium value in a way that depends on the algorithm. We define the escape time t esc , an integer, as the time t at which d(t) has increased by ten percent: with γ = 0.1. Our conclusions are robust with respect to the value of γ.

Escape times-scaling theory
The local Metropolis algorithm and the straight ECMC both have an intrinsic parameter, namely the range δ of the displacement set or the chain time τ chain . These two parameters play a similar role. Two limiting cases can be analyzed. For the Metropolis algorithm at small δ, a trajectory spanning a constant distance is required to escape from an εrelaxed Böröczky configuration. As the dynamics is diffusive, we have const = δ √ t esc . For the straight ECMC with small chain times τ chain , the effective dynamics (after subtraction of the uniform displacement), is again diffusive. This leads to: (for small δ, τ chain ).
On the other hand, even for large δ or τ chain , the Markov chain must make a certain number of moves on a length scale ε in order to escape from the εrelaxed Böröczky configuration. In the Metropolis algorithm, the probability for a trial on this scale is ε/δ for the cross-shaped displacement set, and ε 2 /δ 2 for the square-shaped displacement set. For the straight ECMC with large τ chain , all displacements beyond a time ∼ ε (or, possibly, ∼ N ε) effectively cancel each other, because the constraint graph is rigid. This leads to: (for large δ, τ chain ).
The two asymptotes of eqs (10) and (11) form a "V " with a base δ V at ∼ 3 √ ε (for the Metropolis algorithm with a cross-shaped displacement set, and for straight ECMC) and at ∼ √ ε (for the Metropolis algorithm with a squareshaped move set). The resulting optimum, the minimal escape time with respect to ε, is These scalings balance two requirements: to move by a constant distance (which favors large δ or τ chain ) and to move on the scale ε (which favors small δ or τ chain ).
The forward, reflective, and Newtonian ECMC move in any direction, even in the absence of resamplings, so that their displacement sets are 2Ndimensional. This avoids the rigidity problem of straight ECMC (the fact that the number of constraints can be larger than the number of degrees of freedom). These algorithms introduce no intrinsic scale (as δ or τ chain ). The effective step size of moves may thus adapt as the configuration gradually escapes from the ε-relaxed Böröczky configuration. The step size is initially on the scale ε, but then grows on average by a constant factor at each event, reaching a scale ε > ε after a time ∼ ln(ε /ε). The scale ε at which the algorithms break free is independent of the initial scale ε, and we expect a logarithmic scaling of the escape time (measured in events): t esc ∼ ln(1/ε) (reflective, forward, and Newtonian ECMC).

Escape times-computation results
We now test the scaling theory (see Subsection 3.2.2) of the escape times for k = 5 ε-relaxed Böröczky configurations for small relaxation parameter ε. For the local Metropolis algorithm and the straight ECMC, the predicted behavior of t esc for small and large parameters δ or τ chain is clearly visible (see eq. (12) and Fig. 4).
The absence of an imposed scale for displacements manifests itself in the forward ECMC in the logarithmic dependence on time of the mean free path, that is, the ensemble-averaged displacement between events. As the velocity has unit value, the free path is equal to the difference of Monte-Carlo times t MCMC (t + 1) − t MCMC (t) between subsequent events. Individual evolutions as a function of time t for small relaxation parameters ε and ε nicely overlap when shifted by their escape times (see Fig. 5). The time t here refers to the number of events and not to the Monte-Carlo time t MCMC , which depends exponentially on the number of events t. Starting from an ε-relaxed Böröczky configuration with ε = 10 −30 , as an example, the same number of events is on average required to move from a mean free path of ∼ 10 −30 to 10 −25 , as from a mean-free path ∼ 10 −25 to 10 −20 (see Fig. 5). Overall, escape times (with optimized intrinsic parameters for the Metropolis algorithm and for straight ECMC), validate the algebraic scalings of eq. (12), on the one hand, and the logarithmic scaling of eq. (13), on the other (see Fig. 6). Newtonian ECMC appears a priori as the fastest variant of ECMC. However, it frequently gets gridlocked, i.e., trapped in circles of repeatedly active disks with a diverging event rates. Gridlocks also rarely appear in straight and reflective ECMC. In runs that end in gridlock, escape times are very large, possibly diverging (in Figs 4 and 6, median escape times are therefore displayed for these algorithms, rather than the means). The fraction of gridlocking simulations increases with 1/ε. For the Kahle core, this effect is negligible for all ε. For the Böröczky core, Newtonian ECMC runs into gridlock for roughly one third of individual simulations for ε = 10 −29 (see Fig. 6b, the logarithmic scaling is distorted even for the median). For further discussion of gridlocks, see Subsection 4.1.3.

Discussion
In the present section, we discuss our results for the escape times (Subsection 4.1), as well as a number of more fundamental aspects of Böröczky packings in the context of MCMC (Subsection 4.2).

Analysis of measured escape times
ECMC is a continuous-time MCMC method, and its continuous time t MCMC takes the place of the usual count of discrete-time Monte-Carlo trials. In ECMC, each chain corresponds to a segment [t MCMC , t MCMC +τ chain ) of Monte-Carlo time. However, ECMC is event-driven. The time t, and especially the escape time t esc , are integers, and they count events. The computational effort in hard-disk ECMC is O (1) per event, using a cell-occupancy system that is also implemented in the BigBoro software package. In several of our algorithms, the times t and t MCMC are not proportional to each other, because the mean-free path (roughly equivalent to the time between events) evolves during each individual run.

Range of speedups
The speedup realized by lifted Markov chains, of which ECMC is a representative, corresponds to the transition from diffusive to ballistic transport [25,26,6]. For Markov chains in a finite sample space Ω, the Monte-Carlo time for mixing of the lifted Markov chain cannot be smaller than the square root of the mixing time for the original (collapsed) chain. The remarkable power-lawto-logarithm speedup in ε realized by some of the ECMC algorithms concerns times which measure the number of events. The Monte-Carlo escape times probably conform to the mathematical bounds, although it is unclear how to approximate hard-disk MCMC for ε → 0 through a finite Markov chain. Mathematical results for the escape times from locally blocked configurations would be extremely interesting, even for models with a restricted number of disks.

Space of ε-relaxed Böröczky configurations
Any ε-relaxed Böröczky configuration is merely a sample in a space B ε of volume ∼ ε 2N . We have in fact checked that the position x i of disk i in that configuration can be replaced by x i +ε i (whereε i a random vector inside the circle of radius ε) without affecting the scaling of escape times expressed in eqs (12) and (13). For a given upper limit t cpu of CPU time, this corresponds to a volume of B ε (that cannot be escaped from in t cpu ) scaling as ∼ t −3N cpu , for example, for the straight ECMC and scaling as ∼ exp (−2N t cpu ) for the forward ECMC. We expect B ε to have a triple role, as a space of configurations containing bottlenecks (see Subsection 4.1.3), as a space of configurations that the Monte-Carlo dynamics cannot practically escape from, but maybe also as a space that it cannot even access.
The volume of "practically" stable configurations, as well as the corresponding changes in the free energy per disk are probably unmeasurably small. It is however remarkable that these excluded volumes cannot be escaped from. In many MCMC algorithms for physical systems, as for example the Ising mode, parts of sample space are practically excluded because of their low Boltzmann weight, but they do not feature diverging escape times at finite N .

Gridlock of hard-disk ECMC algorithms, resamplings
ECMC algorithms for soft potentials require random numbers at each event. In contrast, the hard-disk ECMC algorithms of Subsection 3.1.2 except of the forward ECMC treat events through deterministic collision rules. At high density, this can make them susceptible to gridlock, in other words to diverging event rates of chains with successive disks in permanent contact. The Monte-Carlo time between events then goes to zero. Gridlock plays no role in large systems at reasonable densities, but it has been discussed in straight ECMC [27].
Gridlock is the very essence of ECMC dynamics that starts from a Böröczky packing, but it also appears as a final state for ε-relaxed Böröczky configurations. In runs from such configurations, we observe gridlock mostly for Newtonian ECMC with the Böröczky core, rendering the analysis of its scaling behavior with the relaxation parameter ε impossible. It also rarely appears in straight and reflective ECMC for the smallest ε. Because of the infinite event rate, gridlock cannot be remedied through resamplings after a finite Monte-Carlo time t MCMC . To overcome gridlock, one can probably introduce event-based randomness to Newtonian ECMC as is done in forward ECMC.

Böröczky packings and local MCMC: fundamental aspects
We now discuss fundamental aspects of Böröczky packings, from the issue of irreducibility to the question of statistical ensembles, the connection with bottlenecks and, finally, to non-local MCMC algorithms.

Irreducibility of local hard-disk MCMC
Strictly speaking, ECMC can be irreducible only if Ω/B is connected, where B is a suitably defined space of locally stable configurations. Packings in B (a space of low dimension) are certainly invariant under any version of the ECMC algorithms, so that they cannot evolve towards other samples in Ω. Connectivity in Ω/B would at least assure that this space can be sampled. In addition it appears necessary to guarantee that a well-behaved initial configuration cannot evolve towards B or even towards an ε-environment around it. These two properties appear not clearly established for finite densities η and for large N . (At small N , counter-examples are easy to construct.) In other models, for example the Ising model of statistical physics, irreducibility can be proven for any N .
These unresolved mathematical questions concerning irreducibility do not shed doubt on the practical usefulness of MCMC for particle systems. First, the concept of local stability is restricted to hard disks and hard spheres (that is, to potentials that are either zero or infinite). The phase diagram of soft-disk models can be continuously connected to the hard-disk case [28]. For soft disks, irreducibility is trivial, but the sampling speed of algorithms remains crucial. Second, in applications, one may change the thermodynamic ensemble. In the N P T ensemble (further discussed in Subsection 4.2.2), the central simulation box fluctuates in size and can become arbitrarily large. In this ensemble, irreducibility follows from the fact that large enough simulation boxes are free of steric constraints. Again, the question of mixing and correlation time scales is primordial. Third, practical simulations that require some degree of irreducibility are always performed under conditions where the simulation box houses a number of effectively independent copies of the system. This excludes the crystalline or solid phases. Monte Carlo simulations of such phases are more empirical in nature. They require a careful choice of initial states, and are then not expected to visit the entire sample space during their time evolution. Fundamental quantitative results can nevertheless be obtained [29].

Böröczky packings and the N P T ensemble
The concepts of packings and of local and collective stability make sense only in the N V T ensemble, that is, for a constant number of particles and for a simulation box with fixed shape and volume (the temperature T = 1/β that appears in N V T plays no role in hard-disk systems [12]). In the N P T ensemble, the pressure P is constant, and the size of the simulation box may vary. The equivalence of the two ensembles is proven [30] for large N , so that the choice of ensemble is more a question of convenience than of necessity. As we will see, in the N P T ensemble, tiny relaxation parameters (as ε = 10 −29 in Fig. 6) are not maintained for normal pressures and system sizes.
To change the volume at constant pressure, one may, among others, proceed to "rift volume changes" (see [31,Sect. VI]) or else to homothetic transformations of the central simulation box. We discuss this second approach (see [12,Sect. 2.3.4]), where the disk positions (but not the radii) are rescaled by the box size L as: Each configuration is then specified by an α vector in the 2N -dimensional periodic unit square and an associated volume V = L 2 , which must satisfy V ≥ V cut (α). A classic MCMC algorithm [32] directly samples the volume at fixed α from a gamma distribution above V cut (α), below which (α, V ) ceases to represent a valid hard-disk configuration [12, eq. 2.19]. Typical sample volumes are characterized by βP (V − V cut ) ∼ 1, and with V = (L cut + ∆L) 2 , it follows that This equation illustrates that a packing, with ε → 0, is realized as a typical configuration only in the limit βP → ∞. For the Böröczky packings of Fig. 1, we have L 20, and a typical value for the pressure for hard-disk systems is βP ∼ 1, which results in ε ∼ 10 −3 . In the N P T ensemble, as a consequence, escape times from a packing naturally correspond to a relaxation parameter ε ∼ 1/(βP V ), in our example to t esc (ε ∼ 10 −3 ), which is O (1). The above N P T algorithm combines constant-volume N V T -type moves of α with the mentioned direct-sampling moves of V at fixed α. In practice, however, N P T calculations are rarely performed in hard-disk systems [33,34]. This is because, as discussed in eq. (15), the expected single-move displacement in volume at fixed α is ∆V ∼ 1/(βP ), so that ∆V /V ∼ 1/N (because N ∼ V and βP ∼ 1). The fluctuations of the equilibrium volume V eq (averaged over α) scale as √ V eq , which implies ∆V eq /V eq ∼ 1/ √ N . The volume-sampling algorithm requires ∼ N single updates of the volume to go from the 1/N scale of volume fluctuations at fixed α to the 1/ √ N scale of the fluctuations of V eq at equilibrium. This multiplies with the number of steps to decorrelate at a given volume. In practice, it has proven more successful to perform single N V T simulations, but to restrict them to physical parameters where the central simulation box houses a finite number of effectively independent systems mimicking constant-pressure configurations.

Bottlenecks in MCMC algorithms
Markov chains can be interpreted in terms of a single bottleneck partitioning the sample space into two pieces [24,Sect. 7.2]. The algorithmic equilibrium flow across the bottleneck sets the conductance of an algorithm, which again bounds mixing and correlation times. Ideally, MCMC algorithms would be benchmarked through their conductances.
In the hard-disk model, the bottleneck has not been identified, so that the benchmarking and the analysis of MCMC algorithms must rely on empirical criteria. However, Böröczky packings and the related ε-relaxed Böröczky configurations may well exemplify possible bottlenecks and the escape times studied in Subsection 3.2 may model mixing times. They certainly provide lower bounds. Most importantly, the benchmarks obtained by comparing escape times may carry important lessons on the relative merits of sampling algorithms.

Böröczky packings and non-local MCMC
In this work, we concentrate on local MCMC algorithms, with infinitesimal displacements (for ECMC) or very small displacements (for the Metropolis algorithm), because real-life continuous-space problems usually require the use of local methods [6]. Global-move algorithms, as the cluster algorithms in spin systems, rely on a priori probabilities for many-particle moves that are too complicated. On the other hand, global single-particle moves are related to the single-particle insertion probabilities, in other words to fugacities (the exponential of the negative chemical potential) that are prohibitively small.
In view of the scarcity of exact results for hard-disk MCMC algorithms, we now discuss the global-move Metropolis algorithm in which at each time step a randomly chosen disk is placed at a random position inside the box. This corresponds to the Metropolis algorithm of Subsection 3.1.1 with a squareshaped displacement set and a range δ = L/2. This non-local algorithm has no problem escaping from a Böröczky packing. Moreover, it is proven to mix in O (N log N ) steps at densities η < 1/6 [35,36] (see also [37]). This result has made it possible to prove that the liquid phase in the hard-disk system extends at least to the density η = 1/6 [36]. The density bound for the algorithm (which yields a bound for the stability of the liquid phase) is much smaller than the empirical density bound for the liquid phase, at η 0.70. At this higher density, the global-move Metropolis algorithm and the more general hard-disk cluster algorithm [38] are almost totally stuck. For applications, we imagine Böröczky packings to be part of configurations at such high densities, where global moves cannot be used.

Conclusion
Building on an early breakthrough by Böröczky, we have studied in this work locally stable hard-disk packings. Böröczky packings are sparse, with arbitrarily small densities for large numbers N of disks. We constructed different types of these packings to arbitrary precision for finite N and made our implementation openly accessible. Böröczky packings are locally, but not collectively stable. Using singular-value decomposition (in an implementation that is included in our open-source software) we explicitly exposed the unstable collective modes. We furthermore reduced the radius of Böröczky packings slightly, and determined the escape times from ε-relaxed Böröczky configurations as a function of the parameter ε for a number of local MCMC algorithms, including several variants of ECMC, arbitrary-precision implementations of which are also made openly available. Although the algorithms depart from each other in seemingly insignificant details only, we witnessed widely different escape times, ranging from 1/ε to log(1/ε). Our theory suggested that the significant speedup of some of the algorithms is rooted in their event-driven nature coupled to their lack of an intrinsic scale. We pointed to the importance of statistical ensembles to reconcile the obvious loss of irreducibility in the presence of Böröczky packings with the proven practical usefulness of local hard-disk MCMC algorithms.
We expect the observed differences in escape times to carry over to realworld ECMC implementations. In statistical mechanics, bottlenecks and escape times possibly play an important role in polymer physics and complex molecular systems and some of the algorithms studied here may find useful applications. Escape times may also play an important role in the study of glasses and in granular matter, where the high or even infinite pressures favor local configurations that resemble the mutually blocked disks in the ε-relaxed Böröczky configurations. We finally point out that the very concept of locally blocked packings naturally extends to higher dimensions. layers (see Section 2.1.3). The geometric convex polygonal chain A geo with different attenuation parameters φ, and the circular A circ are implemented (see Section 2.2.1). The core, the number of layers, and the convex polygonal chain are specified using command-line arguments that are described in the README.md file of the package, as well as in the output of the script's --help command-line option. The docstrings of the script contain further information.
The Python script construct packing.py uses arbitrary-precision decimal floating-point arithmetic (using the decimal module of Python's standard library). Two additional command-line options specify the number of decimal digits, and the precision b of the bisection search for the value g < 2 that renders the Böröczky packing compatible with periodic boundary conditions (see Section 2.1.3). In general, b should be smaller than the number of places of the decimals. The script modifies If the bisection search succeeds, the script first tests that no pair of disks has a distance smaller than 2 − 10 −b+2 . Second, it checks that every disk has at least three contacts with distances in the interval [2 − 10 −b+2 , 2 + 10 −b+2 ], and, finally, that the total number of contacts agrees with eq. (4). The final configuration and its parameters (as for example the system length) are stored in a human-readable format in a specified output file.
The example packings directory of BigBoro contains several Böröczky packings. The packings are contained in corresponding subdirectories (as for example kahle geometric 5). The headers of these files contain the values of the command-line arguments for construct packing.py. A plot of each example configuration is provided. The different packings in kahle geometric 5 and boro geometric 5 (see Fig. 1) were heavily used in this work. Although the bisection search of the Böröczky-packing construction usually requires an increased precision, the high-precision packings with small enough k may be used as input for standard double-precision applications. For simplicity and improved readability, we provide packing double.txt files that store the configurations with double precision, where applicable.
A.2 Python script collective escape modes.py The double-precision Python script collective escape modes.py identifies the orthonormal basis vectors ∆ a of the escape matrix M esc from a packing x (see eq. (6)) that have zero singular values. This is the solution space for 2Ndimensional displacements ∆ = {∆ x 1 , ∆ y 1 , ∆ x 2 , ∆ y 2 , . . . } that have a vanishing first-order term in eq. (5) and thus for collective infinitesimal displacements ∆ of all disks that escape from the packing. The script asserts that the configurations x + 10 −8 ∆ a are without overlaps and that all contacts persist at a precision 10 −8 (contacts are lost at second order only). Furthermore, the script ensures that the uniform translations of all disks along the x-and yaxis are part of the solution space. Finally, the basis vectors ∆ a are stored in a human-readable output file, and optionally represented as in Fig. 2. The input filename of the packing, and the output filename for the collective escape modes are specified in command-line arguments. Further optional arguments specify the filename for the plots of the escape modes, and the system length of the central simulation box (that is unnecessary for packings generated by the Python script construct packing.py in which case the system length is parsed from the packing file). The package's README.md file, as well as the --help command-line option and the docstrings of the script contain more detailed information.

A.3 Go application go-hard-disks
The Go application go-hard-disks relies on a cell-occupancy system for the efficient simulations of large-N hard-disk systems using several variants of the ECMC algorithm. Straight, reflective, forward, and Newtonian ECMC are implemented. In its current form, it samples the maximum nearest-neighbor distance d(t) (see eq. (8)) after a given sampling time. All computations use a fixed number of mantissa bits (in base 2) that may exceed the usual 24 or 53 bits for single-or double-precision floating-point values. (We use the math/big package of the Go standard library for the arbitrary-precision arithmetic.) The ECMC variant, its parameters (as for example the sampling time or chain time), and further specifications (the number of mantissa bits, the cell specifications, the filename for the initial configuration, etc.) are again set using command-line arguments (see the README.md file of the package for details on the installation process and the possible arguments).

A.4 License, dependencies, software versions
The BigBoro software package is published as an open-source project under the GNU GPLv3 license. It is available on GitHub as part of the JeLLyFysh organization. 1 Users can clone or fork the repository to study the code, and to run the Python3 scripts construct packing.py and collective escape modes.py, and the Go application go-hard-disks. The Python3 scripts rely on NumPy as their only external dependency [39]. Optional plotting also requires the Matplotlib library [40]. They are expected to work with any Python3 version and any NumPy version ≥ 1.20. Python scripts were tested with Python 3.9 and NumPy 1.21. The Go application only requires the Go standard library. It is expected to work with any Go version ≥ 1.13 (tested with 1.16). Users can communicate with the authors (for suggestions or bug reports, etc.) through GitHub issues, and are encouraged to contribute to the project by pull requests.