Voronoi diagrams in quasi-2D hard sphere systems

Variants of the Voronoi construction, commonly applied to divide space, are analysed for quasi-two-dimensional hard sphere systems. Configurations are constructed from a polydisperse lognormal distribution of sphere radii, mimicking recent experimental investigations. In addition, experimental conditions are replicated where spheres lie on a surface such that their respective centres do not occupy a single plane. Significantly, we demonstrate that using an unweighted (no dependence on sphere size) two-dimensional Voronoi construction (in which the sphere centres are simply projected onto a single plane) is topologically equivalent to taking the lowest horizontal section through a three-dimensional construction in which the division of space is weighted in terms of sphere size. The problem is then generalised by considering the tessellations formed from horizontal sections through the three-dimensional construction at arbitrary cut height above the basal plane. This further suggests a definition of the commonly-applied packing fraction which avoids the counter-intuitive possibility of it becoming greater than unity. Key network and Voronoi cell properties (the fraction of six-membered rings, assortativity and cell height) and are analysed as a function of separation from the basal plane and the limits discussed. Finally, practical conclusions are drawn of direct relevance to on-going experimental investigations.

Abstract. Variants of the Voronoi construction, commonly applied to divide space, are analysed for quasi-two-dimensional hard sphere systems. Configurations are constructed from a polydisperse lognormal distribution of sphere radii, mimicking recent experimental investigations. In addition, experimental conditions are replicated where spheres lie on a surface such that their respective centres do not occupy a single plane. Significantly, we demonstrate that using an unweighted (no dependence on sphere size) two-dimensional Voronoi construction (in which the sphere centres are simply projected onto a single plane) is topologically equivalent to taking the lowest horizontal section through a threedimensional construction in which the division of space is weighted in terms of sphere size. The problem is then generalised by considering the tessellations formed from horizontal sections through the three-dimensional construction at arbitrary cut height above the basal plane. This further suggests a definition of the commonly-applied packing fraction which avoids the counter-intuitive possibility of it becoming greater than unity. Key network and Voronoi cell properties (the fraction of six-membered rings, assortativity and cell height) and are analysed as a function of separation from the basal plane and the limits discussed. Finally, practical conclusions are drawn of direct relevance to on-going experimental investigations.

Introduction
Hard particle models are a central tenet of statistical physics, forming the basis of fundamental research dating from the earliest computations to current research [1]. Despite their simplicity, the hard disk and hard sphere models are able to explain many of the behaviours of classical particles in two and three dimensions. In particular, such models are highly applicable to the study of colloids [2]. As part of the characterisation of these systems, it is useful to divide the sample space into regions associated with each particle. This allows determination of the coordination environments around each particle, so that network properties such as the neighbour degree distribution and correlations can be calculated [3][4][5][6]. This information can be used for example to give insights as to the phase behaviour in these systems [7][8][9][10].
The division of space into appropriate cells is most commonly achieved by construction of a Voronoi diagram [11]. In the simplest unweighted version, the Voronoi diagram partitions the sample into a system of tessellating polyhedra (or polygons in 2D), which encapsulate all the space closest to each sphere (or disk). The elegance of the Voronoi diagram arises from the fact that only the sphere centres are required for its construction, with no requirement for a cut-off parameter. While the unweighted Voronoi tessellation is very effective for studying monodisperse spheres, there are some limitations for polydisperse species. Specifically, the Voronoi partition underestimates the space assigned to large spheres and overestimates that for small spheres-a simple reflection of the lack of information on the radii. To rectify this, weighted modifications have since been suggested which take account of the differences in radii [12]. With these developments, the Voronoi diagram has proved invaluable in the analysis of theoretical and experimental particle configurations in both two and three dimensions.
However, the recent shift in focus towards development of reduced-dimensional materials has led to investigations into the novel world of quasi-2D hard particle systems. Quasi-2D systems consist of a monolayer of hard spheres either confined between two planes [13,14], adsorbed at an interface [15,16] or sedimented on a surface [17,18]. These quasi-2D arrangements raise an interesting question as how best to perform a Voronoi analysis i.e. whether we can treat the spheres as a simple 2D array of points, or whether it is necessary to calculate a full 3D weighted construction, requiring the radii.
In fact, the nature of the system will determine which Voronoi analysis is most appropriate. If we consider the simplest quasi-2D case, which is a monolayer of monodisperse spheres, the particle centres will occupy the same plane. As demonstrated in figure 1(a), the 3D Voronoi polyhedra for such system are prismatic, such that any horizontal cut through the construction will yield the same 2D tessellation. In addition, it is relatively straightforward to see that this 2D tessellation is topologically equivalent to the 2D Voronoi diagram, calculated in the plane of the spherical particle centres. This is to say that the 3D Voronoi can be trivially constructed from the 2D Voronoi, and so such systems can be analysed with an unweighted 2D Voronoi, with no loss of information.
A simple extension is then to take polydisperse spheres in which the particle centres occupy the same plane (e.g. particles suspended at an interface). As previously alluded to, a weighted Voronoi construction is required in this case, to avoid the possibility of the polyhedra cutting through the larger spheres. Calculating the 3D Voronoi diagram, weighted by the sphere radii, once again leads to prismatic polyhedra, as shown in figure 1(b). Now the 2D tessellations formed from horizontal cuts are topologically equivalent to the 2D Voronoi diagram, weighted by the sphere radii, calculated in the plane of the sphere centres. Once again the partitioning of space can be fully described in two dimensions, and so such systems can be analysed with a weighted 2D Voronoi.
The more complex case is to consider a monolayer of polydisperse spheres sedimented on a surface. Under these conditions the spheres share a basal tangent plane, and the centres no longer occupy a common plane. Although there is a trivial projection of the particle centres into 2D, the interaction distances between spheres are non-additive [19], reflecting the essential 3D nature of the problem. As can be seen in figure 1(c), the polyhedra in the 3D Voronoi diagram, weighted by sphere radii, are no longer prismatic. The 2D tessellations formed from horizontal cuts through the system therefore have a non-trivial relationship with cut height. Rather, these 2D tessellations show an evolution in structure as the height above the basal plane increases-both in terms of the number of polygons and their properties. Under these quasi-2D conditions, it is not necessarily clear what is the best approach for conducting a Voronoi analysis, or even how to define even basic properties such as the packing fraction.
It is these final quasi-2D systems which are the focus of this work. The structure of this paper is then as follows. We will first demonstrate that the basal tessellation in the weighted 3D Voronoi is topologically equivalent to an unweighted 2D Voronoi diagram. We will then go on to explore the more general stereology problem, showing that tessellation at an arbitrary cut height above the basal 2D plane can be related to a 2D weighted Voronoi diagram. These tessellations will also be compared to those formed from equivalent arrangements of hard disks. Finally these conclusions will be tested using numerical Monte Carlo simulations, and conclusions drawn for use in an experimental context.

Quasi-2D hard sphere systems
In this work we consider a system of polydisperse hard spheres sedimented onto a surface, such that all spheres share a common basal tangent plane. For numerical simulations, we will further assume that the radii, R i , of the spheres follow a lognormal distribution: where as usual μ and σ are respectively the mean and standard deviation of the logarithm of the radii. This distribution is chosen to ensure the radii of randomly generated spheres are always positive. While this does not affect the fundamental conclusions of the Voronoi analysis, it will help quantify properties of the system such as the packing fraction.

Weighted Voronoi tessellations
A Voronoi diagram is constructed through the placement of dividing planes, between the centres of neighbouring particles. The intersection of these planes forms the characteristic tessellating polyhedra, each bounding the space associated with a given particle.
In the common unweighted approach, the dividing plane between two neighbouring particles separated by the Euclidean distance d ij , is simply located midway between the particles at a distance d ij /2. To construct a weighted Voronoi diagram, one simply adjusts the position of the dividing plane, such that it is further from the particle with the greater weight. The weighting method we will be using is the so called radical tessellation, or power diagram [20,21]. In this modification, the dividing plane between . Panel (a) shows a monodisperse system with spheres sharing a basal plane and their centres also occupying the same horizontal plane. An unweighted 3D Voronoi diagram is overlaid. The 3D Voronoi poyhedra are prismatic such that any horizontal cut is identical and topologically equivalent to a 2D unweighed Voronoi, calculated in the plane of the sphere centres. Panel (b) shows a polydisperse system where sphere centres lie in the same horizontal plane. A 3D Voronoi diagram weighted by the sphere radii is overlaid. Again the 3D Voronoi poyhedra are prismatic such that any horizontal cut is identical and topologically equivalent to a 2D Voronoi weighted by the sphere radii, calculated in the plane of the sphere centres. Panel (c) shows a polydisperse system where spheres share a basal plane. A 3D Voronoi diagram weighted by the sphere radii is overlaid. Here the tessellations formed from horizontal cuts evolve with cut height, z. The tessellation on the basal plane is topologically equivalent to the 2D unweighted Voronoi diagram, calculated using the projected sphere centres. two particles is placed at a Euclidean distance, d i , from the centre of particle i, given by: where w i and w j are the weights for each particle. The benefit of this method is that is adjusts the partitioning of space so that greater volume is assigned to the particles with larger weight, and is well designed so that all of the sample space remains accounted for-unlike some alternative constructions [22]. In terms of the particle weights, the logical choice in the 3D radical tessellation for a collection of spheres is simply the sphere radii. This is because at the contact distance, d ij = R i + R j , we find from equation (2) that d i = R i i.e. the radical plane sits between the two spheres, producing the most equitable distribution of volume (see for example spheres A, C in figure 2(b)). Furthermore, when the radii are equal, we find d i = d ij /2 and the result from the standard unweighted Voronoi is regenerated as expected.

Stereological relationships
When studying quasi-2D colloids, it is possible to construct a 3D weighted Voronoi using the coordinates of the sphere centres, c i = (x i , y i , R i ). However, in most cases researchers prefer to analyse configurations using a 2D Voronoi using the projected sphere centres c i = (x i , y i ). The rational behind this is clear: 2D analysis is easier . Panel (a) shows two spheres and the dividing radical plane between them. The radical plane intersects the tangent plane at half the horizontal distance between the spheres. This intersection generates the same dividing line in the tangent plane as would the unweighted Voronoi using the projected sphere positions. Panel (b) shows three spheres and the 3D weighted Voronoi polyhedron formed around sphere A. The height of the cell is given by the topmost point, denoted h S (R A ). Panel (c) shows the same three spheres and the polyhedron that would be formed from stacking 2D Voronoi tessellations with disk radii as weights. As the weight for A is not defined above the sphere diameter, the polyhedron is truncated in comparison with in panel (b). The height of the cell is then given by to visualise and rationalise. In addition, as in figure 1(a), for monodisperse systems the 3D representation contains no more information on neighbouring interactions and assigned volumes than the 2D analogue. However, as in figure 1(c), for polydisperse systems while this may be a good first approximation, this is not necessarily accurate, and so it is important to examine the relationship between Voronoi diagrams in two and three dimensions to ensure correct physical meaning is attributed to the Voronoi analysis.

Limiting case with basal plane.
We first consider the stereological relationship between the 3D weighted Voronoi and the 2D tessellation formed from the intersection of the radical planes with the basal tangent plane. In figure 1(c), this corresponds to the bottommost 2D tessellation. To do this take the arrangement in figure 2(a), with two spheres A, B separated by the dividing radical plane, V , and sharing the tangent plane, T , with equation z = 0. The distance between the sphere centres is d AB , while the distance between the projected centres is denoted d AB . The dividing plane V , has a normal given by n = Defining the direction of x as the projection of d An on the basal plane, the plane equation for V can be deduced as: The intersection of V with T occurs at z = 0, and so we find that line of intersection is given simply by x = d AB /2. However we notice something interesting about this result. This dividing line in 2D is the same as that which would be obtained from constructing the 2D unweighted Voronoi using the projected sphere centres c i . As this must be true for all pairs of neighbouring spheres, it follows that the entire 2D unweighted Voronoi will be constructed. This leads to the key relationship of this work: for quasi-2D systems of spheres, the unweighted 2D Voronoi diagram is topologically equivalent to the tessellation formed by taking the basal faces of the polyhedra in the weighted 3D Voronoi diagram. This holds when the 2D Voronoi is constructed using the point set of the sphere centres projected onto the basal plane, and the 3D Voronoi using the point set of the original sphere centres and weighted by the sphere radii.

General case with arbitrary horizontal plane.
Taking the basal faces of the 3D Voronoi polyhedra can equally be thought of as taking a horizontal cut through the tessellation at z = 0. The result above is essentially a special case of the fact that a cut through a 3D Voronoi must yield a tessellation which is equivalent to some weighted 2D Voronoi [23,24]. Following from the result above, we may ask what is the analogous relationship between a 2D and 3D Voronoi when taking a cut at arbitrary z, as in figure 1(c). Revisiting the simple arrangement in figure 2(a), for a given horizontal cut height, z, we now define the projected sphere centres as c i = (x i , y i , z). Regardless of the value of z, we note that the distances between the projected centres remain the same at d AB . To obtain the same distance between the projected sphere centres and the dividing plane using both a 2D and 3D Voronoi, it can be seen from combining equations (2) and (3) that: The simplest solution for the weights that satisfy this equation is therefore: Again as this will naturally extend to a collection of many spheres, we obtain a more general relationship: for quasi-2D systems of spheres, the tessellation formed from a horizontal cut through a 3D Voronoi diagram weighted by the sphere radii, is topologically equivalent to a 2D Voronoi diagram weighted according to equation (5). This holds when the 3D Voronoi is constructed using the point set of the original sphere centres, and the 2D Voronoi using the point set of the sphere centres projected onto the horozontal cut plane. The caveat to this result is that if the cut is height is greater than the maximum cell height for a given sphere, such a sphere must be excluded from the 2D Voronoi calculation. We see that the key result in the above section is simply a limiting case of this more general result, when z = 0 and the weights are trivially zero. In addition, in the limiting case, it follows that all spheres must have a cell in the 2D Voronoi. We shall refer to this 2D weighted Voronoi diagram as the sphere-weighted Voronoi and denote variables associated with it with a subscript S.

Connection to system of hard disks.
A horizontal cut through a 3D system of hard spheres will produce a 2D system of hard disks. The radii of these disks will be related to the radii of the spheres, and if they are used as weights in a 2D Voronoi analysis they will also satisfy equation (4): This suggests that the tessellation formed from the horizontal cut through the 3D Voronoi should also be equivalent to the 2D Voronoi using these disk radii. The same configuration is used as in figure 1, with the horizontal cuts corresponding to the left hand diagrams. Black circles indicate the particle section radii at the cut height, which also correspond to the weights in the disk-weighted Voronoi. Black points indicate the projected particle centres, while grey points indicate particles without corresponding cells in the 2D tessellation.
However, there is a caveat to this result. It is clear that equation (6) is only well defined for z 2R i i.e. the cut height is less than the sphere diameter. There is no such constraint for the polyhedral cells in the 3D weighted Voronoi, which may readily extend above the associated particle. Therefore, the result above is only strictly true when the cut height is less than the smallest sphere diameter. A modification is made where points which have diameters less than the cut height are excluded from the 2D calculation. These excluded spheres will therefore not have cells in the Voronoi tessellation. In this case some 3D cells will become truncated in the 2D representation, as shown in figures 2(b) and (c). This will lead to an increasing difference between the two types of partitions as the cut height is increased, as demonstrated in figure 3. This 2D weighted diagram we shall refer to as the disk-weighted Voronoi and denote variables associated with it with a subscript D.

Properties of Voronoi tessellations
We will discuss important properties of the sphere and disk-weighted Voronoi diagrams which will aid with their analyses. These include the expected number of points at a given cut height, the weight distribution of the included points and the nearest-neighbour distances between them. We note that these are the quantities that appear in equation (2) and therefore govern the wider properties of the Voronoi diagram. We will show the analytic results for our model of lognormally distributed sphere radii, but the results extend to any configuration where the sphere radii are known.

Average nearest-neighbour distances.
A useful metric is the expected distance between points in a Voronoi diagram at different cut heights, as this can be used to rationalise much of the observed behaviour. In order to do this we will assume that spheres are uniformly distributed (the ideal gas), neglecting the short range liquid structure. While this might seem an over-simplification, not only does the approximation get better at lower packing fractions but also as the cut height is increased, spheres are lost from the tessellation and their positions become effectively decorrelated.
First we wish to find the average distance between all points in neighbouring Voronoi cells at a given cut height, d z . We introduce the number density at a given cut height, denoted ρ z . This will be the number of points at height z, divided by the total area of the 2D Voronoi diagram. It can therefore be expressed: where ρ 0 is the number density considering all points (as must be the case at z = 0), and N (z) is the proportion of spheres included at a given cut height. By dividing the area equally between all particles, we expect the average distance between neighbouring points in the 2D Voronoi diagram to follow: Alternatively, we can opt to find the average distance to the nth nearest neighbour, d n z . If again a uniform distribution of spheres is assumed (i.e. the dilute limit), equation (12) of Bhattacharyya and Chakrabarti [25], shows that in two dimensions, this is given by the equation:

Properties of disk-weighted Voronoi.
The properties of the disk-weighted Voronoi are somewhat easier to see and so we begin with these. The height below which a Voronoi cell is defined for a sphere of a certain radius is given exactly by h D (R) = 2R, as in figure 2(c). Hence, for a given cut height, the spheres with associated cells in the tessellation will be all those with radii in excess of the inverse function h −1 D (z) = z/2. The proportion of particles with associated cells above a given cut height, denoted N D (z), is therefore simply related to the cumulative distribution function of the sphere radii: as all spheres with R < z/2 are neglected. The moments of the weight distribution function at a given cut height can then be calculated by integrating over the range of the remaining sphere radii, as required. For example, the analytic solution for the second moment can be found in section 2.5.

Properties of sphere-weighted Voronoi.
For the sphere-weighted Voronoi, the maximum height of the Voronoi cell for a sphere of a certain radius, h S (R), is not precisely defined with R. We will derive an approximate form below, but first note that regardless of the functional form of h S (R), in a similar manner to above we can write an equation for the proportion of spheres with associated cells above a given cut height, denoted N S (z), as: with the moments of the weights at a given cut height provided by: in analogy with the disk-weighted case. We now discuss the functional form of h S (R). Rather than an exact relationship, we wish to find an approximate form for the expected Voronoi cell height for a given sphere radius. The height of a Voronoi cell is a complex function that must depend on both the radius of the associated particle and the distances to neighbouring particles. In the following arguments it will be useful to refer to figures 1-3. Consider a spherical reference particle of a given radius, R . Crucially, the dividing planes forming the Voronoi polyhedron around it will only converge to a point when the neighbouring spheres are larger than the reference particle. For the purposes of this analysis, we can therefore think of the particle operating in a reduced density system containing only the particles with radii greater of equal to the reference particle. The density is therefore,

J. Stat. Mech. (2020) 093201
Voronoi diagrams in quasi-2D hard sphere systems and mean average sphere radius, In addition, the average nth nearest neighbour distance will be given by equation (9), substituting the density from equation (14) above. We therefore consider the average environment around the reference particle to consist of spheres with radii R at located at successive distances d n . These spheres will generate dividing Voronoi planes according to equation (3). If we then assume that on average the highest point in the Voronoi polyhedron will occur above the sphere centre (i.e. x = 0, see for example the cells in figure 1), we find: These can be averaged over m nearest neighbours to give the result: As will be shown in section 3.1 and figure 5, we found that averaging over m = 3 nearest neighbour distances gave a good fit to the observed numerical distribution. We reason this on the basis that the top vertex is formed from the intersection of three neighbouring planes, which originate from the closest three spheres. We note again that this proposed form is not intended to be an exact theoretical result, but rather a possible rationalisation of numerical results which will be discussed later.
The important thing about this functional form is that the expected height of the Voronoi cell scales quickly with sphere radius, unlike in the disk-weighted case. This is a result of the density of larger spheres decreasing quickly with sphere size, such that the average interaction distances become increasingly long and the radical planes approach vertical.

Packing fraction
One complication with quasi-2D hard spheres is that the definition of packing fraction is ambiguous. The 2D packing fraction, φ 2D , is usually defined as where ρ is the number density. For a quasi-2D system this definition is not rigorous, as it can lead to a packing fraction greater than unity. This can be considered an example of the Holmes effect, first noted in the analysis of thin sections of granular media [26,27].
As an example, consider the binary crystalline lattice in figure 4, consisting of a two sphere types with radius ratio γ. At γ = 0 the lattice is the hexagonal close packed structure which has the well known maximal monodisperse packing limit of π 2 √ 3 ≈ 0.907. As γ increases, the smaller spheres swell in the tetrahedral holes without affecting the . Variation in packing fraction for a binary crystal with varying radius ratio, γ. Using the standard 2D packing fraction leads to values greater than unity whereas taking the maximum packing fraction with cut height limits the packing fraction to the known maximal packing limit. large sphere positions, and so the naïve two-dimensional packing fraction rises, peaking at 11π 18 √ 3 ≈ 1.108 when γ = 1 3 (the contact distance). Beyond this radius ratio the small spheres can no longer be accommodated in the tetrahedral holes, and so the larger spheres are forced apart which leads to an initial decrease in packing fraction before the hexagonal close packed limit is approached once more. We see that over this range the packing fraction exceeds unity and therefore obscures the physical meaning.
Moving to 3D does not necessarily solve the issue, as it is now not clear what volume (and therefore number density) is most appropriate. However, considering quasi-2D packings as a series of stacked sections through the 3D system allows an alternative definition of packing fraction: as a function of z, the horizontal cut height, where ρ z and w 2 D z are provided by equations (7) and (11) respectively. The rationale for this equation is as follows. As the value, w 2 D z , measures the average square radius of the circular sections in a cut of given height; here φ (z) quantifies the proportion of the total area occupied by the particle sections in the cut plane.
We propose that the overall packing fraction could thereby by quantified in these systems by max [φ (z)]. Referring again to figure 4, we see that this definition gives values for the packing fraction which are now consistent with intuition. For the region where γ 1 3 and the large sphere lattice is unperturbed by the small spheres, the maximal monodisperse packing limit is found. Again beyond this the packing fraction then drops as the smaller spheres swell above the contact distance before the hexagonal close packed lattice is once again recovered. In addition, if this definition were applied to the simpler systems in figures 1(a) and (b), it would match with packing fraction calculated using the disks of the same radii as the spheres. For a polydisperse system, the explicit form of the second moment for the lognormal distribution can be calculated as,

J. Stat. Mech. (2020) 093201
Voronoi diagrams in quasi-2D hard sphere systems We note the pleasing similarity between this result and equation (6), and see that there will be parity when σ = 0 i.e. the monodisperse limit. Intuitively, the form of the packing fraction with z and its maximal value is independent of the number density and is instead solely a property of the sphere radii distribution.

Network properties
To evaluate any differences between the various tessellations in this work, one can evaluate the network properties [6,28]. A network is formed by assigning a node to each Voronoi cell, and linking the nodes of adjacent cells. The node degree, k, is then equivalent to the number of neighbours of each cell or equally the number of cell edges.
The mean node degree must be equal to six, k = 6, as a result of Euler's formula for connected planar graphs. However, the node degree distribution, p k , (also known as the ring statistics) is able to vary providing the above constraint is satisfied. Hence the proportion of any given node degree can provide an insight into the structure of a tessellation. This is usually chosen to be the proportion of hexagons, p 6 , as this is typically the most prevalent ring size and has a physical reference as p 6 = 1 is the crystalline limit.
A further measure is provided by the assortativity [29]. The assortativity describes the tendency of low degree nodes to be adjacent to high degree nodes in a network-which has previously been quantified for example through the Aboav-Weaire law [30]. The assortativity can be measured through the Pearson correlation coefficient: where e jk is the edge joint degree distribution i.e. the probability of nodes of degrees j, k sharing a link. As it is a correlation coefficient, −1 r 1, and the assortativity has well defined limits. The majority of physical networks are found to have r < 0, indicating that high degree nodes are more likely to be found adjacent to low degree nodes than in a perfectly random network (where r = 0).

Monte Carlo simulations
The findings in this work are tested by numerical simulation methods. Systems of hard spheres, at five different number densities in the range ρ = 0.00 → 0.20, are initialised with radii drawn from the lognormal distribution in equation (1), with μ = 0.0, σ = 0.3. Configurations can then be generated using standard Metropolis Monte Carlo techniques [31]. In each calculation P spheres were placed in a 2D periodic box, then equilibrated with 10 5 Monte Carlo cycles with each cycle consisting of P random moves. The moves comprise both random displacements and a proportion of swap moves to aid equilibration of the polydisperse spheres [32]. After equilibration, 10 5 Monte Carlo cycles are performed with sampling every 100 cycles. Tessellations for each sample configuration can be made on-the-fly using both a weighted 2D Voronoi and the full 3D weighted Voronoi cut at a given height [33]. The periodicity in the system ensures that the Voronoi diagrams generated have no boundary, removing any potential edge effects. This maintains the mean polygon size of k = 6, and ensures all node degrees are satisfied in the calculation of the assortativity. We note that because there are differing number of spheres in the tessellations at a given cut height, care must be taken to ensure comparable statistical sampling for each state point. As discussed previously, the greater the cut height, the fewer spheres there are with associated cells in the tessellation. Therefore to ensure an equal number of cells are averaged for each value of z, simulations must be repeated with a different number of random seeds, proportional to 1/N (z). This sampling also has an impact on the optimal number of spheres to include. If this number is too low the statistics will be affected at high z, regardless of number of starting seeds (consider the case where only a handful of spheres are included in a 2D Voronoi, it becomes impossible to achieve a polygon with a large number of sides). One must therefore balance these statistical considerations with the limits of computational efficiency. In this work we elected to use P = 1000 → 5000 spheres, depending on cut height so that all tessellations contained at least 50 cells, with a total of ∼10 7 cells sampled for each state point.

Numerical results
Results from numerical Monte Carlo simulations are detailed in this section. We will compare different Voronoi methods: (a) The section of a 3D Voronoi weighted with sphere radii, R i , cut at height z. (b) The sphere-weighted 2D Voronoi weighted according to equation (5). (c) The disk-weighted 2D Voronoi weighted according to equation (6).
As explained in section 2.3 methods 1 and 2 are equivalent and produce identical results. Method 3 on the other hand is expected to approximate to the other methods in the region near z = 0.
We begin by discussing the number of cells that can be found in each tessellation, as the cut height is increased, as shown in figures 5(a)-(c). This is measured by the quantity N (z), the proportion of the total spheres that have associated cells in the Voronoi. For the sphere-weighted Voronoi, figure 5(a), we find that the number of cells at a given cut height scales with the system density, such that a plot of N S (z) against z × ρ 0 produces a near universal curve. This curve is fit well by equation (12) when m = 3 i.e. averaging over three nearest neighbours. The curves for m = 2, 4 are also presented, lying either side of the m = 3 curve, systematically under-and over-estimating N S (z). The main small discrepancy between systems of different density comes at the low cut heights, where the effects of short range ordering are still observable, as in figure 5(c). After the large initial drop in the number of cells at small cut heights, we see a long tail in the function N S (z). This is indicative of the fact that the cells for the largest spheres are very persistent. As these spheres are in low density and likely spaced far apart, the dividing planes are near-vertical, leading to a slow convergence of the cell walls.
For the disk-weighted Voronoi, the behaviour in N D (z) is exactly described by equation (10). In this case, the functional form is dependent only on the underlying sphere distribution, and not the system density. In contrast to the sphere-weighted case, the number of cells featured in these tessellations decays rapidly to zero beyond R , simply reflecting that the maximum cell height is defined by the sphere diameter, and so the polyhedra are truncated instead of extending above the sphere. Figure 5(d) explores the change in packing fraction as a function of the cut height, as explained in section 2.5. As can be seen the fundamental shape is invariant to the number density, with ρ 0 merely acting as a scaling factor. Once again there is excellent agreement between the numerics and the theoretical results in equation (19). Table 1 shows how the maximum packing fraction differs from the naive 2D packing fraction that can be calculated from the sphere radii and the cell area as in equation (18). Calculating packing fraction in this Table 1. Difference between the maximum packing fraction as a function of height above the surface, and the 2D packing fraction using sphere radii for different number densities. way always acts to reduce the overall number, and prevents the packing ever exceeding unity. We now examine how these differences in dividing space affect the network properties of the system. Previous studies which have compared the tessellations formed from cuts through 3D packings of mono-and bi-disperse spheres have found differences when comparing the results of the different methods [34,35]. For our two Voronoi methods, we see that there is good agreement at low cut heights, but fundamental differences in the asymptotic limit as the cut height is increased. We compare two network metrics for each Figure 7. Comparison of the inter-particle distances (panels (a) and (b)) and particle weights (panels (c) and (d)) for the sphere-and disk-weighted Voronoi at increasing cut heights. In the sphere-weighted case both measures scale as z 1/2 in the high cut height limit, whereas in the disk-weighted case the inter-particle distance increases exponentially, making the weighting effect negligible. system: the proportion of hexagons, p 6 , and the assortativity, r, displayed in figure 6. As expected we see good agreement between methods in both metrics at low cut height, with convergence in the limit of z → 0. However, at larger z key differences emerge. The disk-weighted Voronoi approaches the 2D Poisson-Voronoi (PV) limit, p 6 ≈ 0.295 and r ≈ −0.15. This is the limit which is obtained from having unweighted points randomly located in 2D space. In contrast the sphere-weighted Voronoi approaches a different, less well understood, limit with more diverse ring distribution and reduced ring-ring correlations. We will discuss the nature of this limit in a following paragraph.
To explain these observations we examine two further measures, the expected weight at a given cut height, w z , and the average nearest neighbour distance at a given cut height, d z , given in figure 7. We will refer to these in reference to equation (2), which defines the position of the radical plane between two particles. For the disk-weighted Voronoi, we showed above that the system approached the random PV limit. In this limit the particles are randomly positioned and effectively unweighted. The only way this is achievable is if the inter-particle distances scale in excess of the weightings with increasing cut height. As seen from figures 7(b) and (d), this is indeed the case, with the distances increasing exponentially as the cut height exceeds an increasing number of particle diameters while the average weighting remains relatively constant. Two further observations arise from these plots: equation (8) describes the average distances very well (with the greatest deviation for higher density systems reflecting the presence of short range liquid structure), and the form of w D z mirrors the oscillations in the network properties. These oscillations, which are particularly pronounced in p 6 ( figure 6(b)), are a result of the balance of weightings with cut height which become increasingly visible at higher densities where the inter-particle distances are smaller.
Applying the same reasoning to the sphere-weighted Voronoi, it is initially difficult to see how the network properties can tend to a limit which is not PV. However, this type of behaviour is not unknown, being similar to the observation that taking random sections through 3D PV polyhedra leads to 2D polygons which follow a distribution other than PV [36]. We have said that as cut height is increased, particle cells are removed at random from the tessellation, so we must be approaching some random limit. Examining equation (2) again, the only way this cannot be PV is if the particle weightings scale in the same way as the inter-particle distances. Figures 7(a) and (c) show that this is indeed the case, with both metrics scaling as z 1/2 . As such, the decreasing density of points is directly offset by the increase in weighting and a random limit other than PV is reached.

Experimental interpretation
Having examined the numerical results, we consider what practical conclusions may be drawn from this work for the experimentalist. As mentioned previously, for those studying quasi-2D systems experimentally, it is often preferable to analyse configurations using particle positions projected into the x − y plane, using a 2D Voronoi method. Not only is this more practicable, it is also leads to straightforward visualisation and analysis. The most common approach is therefore to use an unweighted 2D Voronoi [18,37]. This is understandable, as it requires no knowledge of the particle radii, which maybe prove difficult to accurately determine. With binary or polydisperse spherical particle systems, it was previously unclear that the unweighted 2D Voronoi analysis was in fact wholly appropriate, given that for true 2D systems (e.g. polydisperse disks) it is known to allocate space disproportionately.
However, in this paper we have shown that contrary to what might have been expected, the use of an unweighted 2D Voronoi analysis, for polydisperse spheres sedimented on a plane, is completely valid and well defined. The unweighted 2D Voronoi is topologically equivalent to the tessellation formed by taking the basal faces of the polyhedra formed from the fully weighted 3D Voronoi diagram. Therefore, in the absence of the information about sphere radii, the unweighted 2D Voronoi can still be mapped to the basal section through the weighted 3D diagram. In addition, for such quasi-2D systems, the unweighted 2D Voronoi diagram may be considered the best construction possible. This is because it corresponds to the only section through the 3D weighted Voronoi which is guaranteed to include a cell for each particle; notwithstanding the fact that it is also the simplest Voronoi method to implement and evaluate. We note as an aside, that this section also seems to maximise the proportion of six membered rings in the tessellation. Furthermore, we have also considered the stereology problem for quasi-2D systems, which is a subject of interest in similar studies of polycrystalline materials [38,39]. Here we have shown that the horizontal 2D sections through the 3D weighted Voronoi diagram can be calculated using a 2D weighted Voronoi, with weightings given in equation (5). Therefore, if the sphere radii are known, one can easily calculate the 2D sections at a given cut height. This is convenient as such computational methods are more readily available than cutting the 3D weighted Voronoi directly, which is a non-standard technique.
Finally we note that an alternative definition of the packing fraction can be adopted by considering the quasi-2D system as a series of hard disk tessellations. This packing fraction is well defined in that it cannot exceed unity and does not require an arbitrary definition of the sample volume.

Conclusion
In this work we have explored the relationships between various tessellations in quasi-2D hard sphere systems of direct significance to on-going experimental investigations. Significantly, we have shown the link between the application of an unweighted twodimensional Voronoi construction and a three-dimensional construction in which the division of space is weighted in terms of particle size, for the case in which the spheres sit on a surface. As a result, we have provided a clear geometrical meaning to the commonly used unweighted 2D Voronoi diagram; showing it to be equivalent to the tessellation formed from taking the basal polygons in the 3D Voronoi diagram weighted by the sphere radii.
In addition we have shown how the sphere packing can be analysed as a function of the separation from the basal plane. This has led to a novel definition of the packing fraction which naturally avoids the fraction becoming greater than one. Furthermore, we have shown how the tesselations formed by taking cuts at different heights above the basal plane are related to key network metrics (here the fraction of six-membered rings and the assortativity). Finally, we have highlighted how the key conclusions drawn are of direct relevance to on-going experimental investigations.