Table Of Content

Large-n conditional facedness m n 4 of 3D Poisson-Voronoi cells 1 0 2 H.J. Hilhorst n a J 0 Laboratoire de Physique Théorique, baˆtiment 210 1 Université Paris-Sud and CNRS, 91405 Orsay Cedex, France ] h January 13, 2014 c e m - t Abstract a t s We consider the three-dimensional Poisson-Voronoi tessellation and . at studytheaverage facedness mn of acell known toneighbor ann-faced m cell. Whereas Aboav’s law states that m = A+Bn−1, theoretical n - arguments indicate an asymptotic expansion m = 8+k n−1/6+.... d n 1 Recent new Monte Carlo data due to Lazar et al., based on a very n o large data set, now clearly rule out Aboav’s law. In this work we de- c termine the numerical value of k and compare the expansion to the [ 1 Monte Carlo data. The calculation of k involves an auxiliary pla- 1 1 nar cellular structure composed of circular arcs, that we will call the v 3 Poisson-Mo¨bius diagram. It is a special case of more general Möbius 6 diagrams (or multiplicatively weighted power diagrams) and is of in- 3 terest for its own sake. We obtain exact results for the total edge 2 . length per unit area, which is a prerequisite for the coefficient k , and 1 1 0 a few other quantities in this diagram. 4 1 Keywords: Poisson-Voronoi diagram, Aboav’s law, Möbius diagram, : v large-n behavior i X r a LPT-Orsay-14-01 1 1 Introduction C ellular structures, or spatial tessellations, are of interest because of theirverywideapplicability. Theperhapssimplestmodelofacellu- lar structure is the Poisson-Voronoi tessellation (or ‘diagram’), obtained by constructing the Voronoi cells around pointlike ‘seeds’ distributed randomly and uniformly in space. Whereas two- and three-dimensional Poisson-Voronoi diagrams are relevant for real-life cellular structures, the higher-dimensional case appears in data analyses of various kinds. An excel- lent overview of the many applications is given in the monograph by Okabe et al. [1]. Beginning with the early work of Meijering [2], much theoretical effort has been spent on finding exact analytic expressions for the basic statistical properties of the Voronoi tessellation, in particular in spatial dimensions d = 2 and d = 3, but also in higher dimensions. Of interest is first of all is the probability p (d) that a cell have exactly n n sides (in dimension d = 2) or n faces (in dimension d = 3). Next comes the conditional sidedness (or facedness), usually denoted m (d), i.e. the average n number of sides (or faces) of a cell known to neighbor an n-sided (or n-faced) cell. There has been considerable theoretical interest in the dependence of p (d) and m (d) on n, but only very few analytic results exist. In this work n n we will be interested in m (2) and m (3). n n In two dimensions experimental data are fairly numerous but usually cover a limited range of n values, not beyond n ≈ 10. The data are most often plotted as nm versus n. In the experimental range it has often been n possible to fit them by what is known as Aboav’s ‘linear’ law [3], which says that nm = An+B, where A and B are adjustable parameters. On the basis n of Monte Carlo simulations [4] it has been known since a long time, however, that two-dimensional Poisson-Voronoi cells violate Aboav’s law, the graph of nm being slightly but definitely curved. n In earlier work [5, 6, 7] we have been interested in Voronoi cells with a very large number n of sides (or faces). We determined the exact asymptotic behavior of p (2) in the large n limit and deduced [8] from it, under very n plausible hypotheses, the asymptotic behavior of m , n 1 m (2) = 4+3(π/n)2 +..., n → ∞, (1.1) n whichrulesoutAboav’slaw. Whentruncatedafterthesecondterm, Eq.(1.1) is in quite reasonable agreement with the Monte Carlo data. An extension [9] of these arguments to higher dimensions, under plausible but unproven assumptions, led to −1 m (3) = 8+k n 6 +..., n → ∞. (1.2) n 1 2 Apart from the precise structure of this formula, its most important predic- tion is that Aboav’s linear law is violated also by three-dimensional Poisson- Voronoi cells. At the time, however, the existing d = 3 Monte Carlo data were insufficiently precise to confirm this. Indeed, three-dimensional Monte Carlo results due to Kumar et al.[10] covering the range 10 ≤ n ≤ 22 were interpreted by Fortes [11] in terms of Aboav’s law. Thesituationhaschanged recently duetoanimpressive largescale Monte Carlo simulation by Lazar et al.[12], which provides a rich trove of infor- mation about the three-dimensional Poisson-Voronoi tessellation. Amidst a wealth of other data the authors determine the values m (3) based on a data n set of 250 million Voronoi cells. Their results clearly show the nonlinearity of nm (3). Given these new data it therefore becomes of interest to consider n againthe asymptotic expansion (1.2) and to try and determine the numerical value of the coefficient k . We do so in this paper and compare the result 1 to the Monte Carlo data of Lazar et al. A juxtaposition of the two- and the three-dimensional m is also illuminating. n Insection2werecallhowthequestionofcalculatingthethree-dimensional m in the large-n limit leads to the problem of a special (non-Voronoi) n tessellation on a spherical surface of radius ∼ n1/3, i.e. essentially a two- dimensional problem. This tessellation, whose edges are circular arcs, is of interest in its own right. It is closely related to the multiplicatively weighted (or: Möbius) diagrams reviewed in Ref.[1], which is why we call it the Poisson-Möbius diagram. Section 3 deals with this auxiliary problem and may be read independently of the rest of the paper. We derive the exact expression for a prerequisite for finding k , viz. the average edge length per unit area in the 1 Poisson-Mo¨bius diagram. In section 4 we briefly describe some Monte Carlo work that we did on this tessellation. In section 5 we the return to the three-dimensional m (3) and provide n extensions of Eq.(1.2). 2 The many-faced 3D Poisson-Voronoi cell We consider a three-dimensional Poisson-Voronoi diagram of seed density ρ. This density may be scaled to unity but we will keep it as a check on dimensional consistency. Let the cell of a central seed have n faces. It was argued in Ref.[9] that in the limit of large n certain cell properties become deterministic, in analogy to what happens in a statistical system in the thermodynamic limit. In particular, in the limit of large n the n first- neighbor seeds F lie in a spherical shell of radius R ≃ (3n/4πρ)1/3 (this j n radius was called 2R∗ in Ref.[9]) and of effective width ∼ n−2/3. For the 3 present purpose this width may be set to zero and for n → ∞ the shell may be approximated locally by a flat plane F as shown in Fig.1. Also in that limit, the Voronoi cells of the first neighbors F approach j prisms that intersect F according to the two-dimensional Voronoi diagram of the set of seeds {F }. There is no reason for these seeds to be Poisson j distributed, but their average sidedness is necessarily exactly six, which is therefore also the average number of lateral faces of a prism. Each prism furthermore has at its lower end a face in common with the central Voronoi cell, notshown inthefigure. Attheir upper ends theprismshave facesincommon withthesecond-neighbor cellsconstructed aroundtheseeds S ,S ,.... Akey 1 2 observationisthatforn → ∞thefirst-neighborseedsbecomeinfinitelydense inF, andthatinthatlimitthesurface(tobecalledΓ)separatingthesecond- neighbor cells from the first-neighbor ones becomes piecewise paraboloidal, the piece P (see Fig.1) lying on the paraboloid of revolution equidistant j from S and from F. The second-neighbor seeds have an n independent j spacing ∼ ρ−1/3 between themselves, whereas the typical diameter of a prism vanishes as ∼ n−1/6. Fig.1showstotheleftthegenericcasewhereafirst-neighborcellarounda seed F has asingle faceatitsupper end. This happens witha probability, to 0 be called f , that tends to unity when n → ∞. The same figure shows to the 8 right the exceptional case where the upper end of a first-neighbor cell around aseedF hastwo facesincommonwiththesecond-neighbor cells. Wedenote 1 the probability for this to happen by f . This event occurs only when the 9 upper end of the prism intersects the joint between two paraboloidal surface segments. In the figure to the right, the arc AB is such a joint, itself located in the plane Q that perpendicularly bisects the vector S −S . In Ref.[9] it 1 2 was argued that f = 1−O(n−1/6) and that f = k n−1/6 +..., whereas the 8 9 1 analogously defined probabilities f and beyond are proportional to higher 10 powers of n−1/6. As a consequence a first-neigbor cell will be, upon averaging over the number of lateral faces, eight-faced with a probability f and nine- 8 faced with a probability f . Fromthe relationm (3) = sf it then follows 9 n s s thatk isalsothecoefficient appearinginEq.(1.2). Atacitandplausible, but 1 P unproven hypothesis, is that cells other than those that are 8− and 9−faced contribute only to higher orede in the n−1/6 expansion that we are about to set up. Our task then is to calculate f to leading order in n−1/6. 9 We now observe that the projection onto the plane F of the set of joints between the P yields a diagram (that we will denote by G and for reasons to i be explained call the Poisson-Möbius diagram), and that the question above amounts to asking which fraction of the Voronoi cells in F is intersected by the edges of G. In section 3 we will mathematically formulate the problem of determin- ing the properties of the projected graph G. We will focus on finding the total edge length per unit area, λ, of G and obtain an exact expression for 4 Q P S P 2 B 2 1 S 1 S 2 P Γ 1 S 1 n−2/3 A F F n−1/6 0 F F1 Figure 1: Both figures show the plane F containing the n ≫ 1 neighbors of a many-faced central cell in a 3D Poisson-Voronoi diagram. The central seed itself is located at a distance R ∼n1/3 below the plane and is not shown. The second- n neighbor seeds, of which S and S are examples, have a density ρ. Left: The 1 2 typical prism shaped first-neighbor Voronoi cell of seed F has its upper end face 0 in contact with a single second neigbor cell. It is therefore eight-faced. Right: The exceptional first-neighbor Voronoi cell of seed F has its upperend in contact with 1 two second-neighbor cells, so that it is nine-faced. Figures taken from Ref.[9]. this quantity. In section 4 we verify our analytic result by a Monte Carlo simulation. Having determined λ we will then in section 5 use it to find a numerical value of k . 1 3 The Poisson-Mo¨bius diagram G 3.1 Definition of G We begin by considering a rectangular box [−L,L]2×[0,L] whose volume we denote by V = 4L3. Let the seeds in this box be located at S ,S ,...,S 1 2 N with N such that N/V = ρ. We will at some convenient point let L → ∞ at fixed ρ, so that the box becomes R3 and the S become Poisson distributed. + i We set S = (x ,y ,z ). i i i i The surface z = P (x,y) given by i z (x−x )2 +(y −y )2 i i i P (x,y) = 1+ (3.1) i 2 z2 (cid:18) i (cid:19) is a paraboloid of revolution of focus S and axis perpendicular to the xy i plane. It separates R3 into a region containing all points closer to the xy + plane than to seed S , and its complement. i Let z = Γ(x,y)be the surface that separates the upper half-space R3 into + a region of points closer to the xy plane than to any of the seeds, and its 5 10 8 6 4 2 0 0 2 4 6 8 10 Figure 2: Poisson-Mo¨bius diagram G obtained by projecting the surface elements P that constitute the surface Γ [see section 3.1] onto the xy plane. j complement. Then Γ(x,y) is built up out of piecewise paraboloidal surface elements P lying on the P . Only paraboloids P associated with seeds S i i i i sufficiently close to the xy plane will contribute surface elements. The arcs along which the surface elements of Γ join will be referred to as ‘joints’. The intersection of two arbitrary paraboloids P and P is an ellipse lo- 1 2 cated in the plane that perpendicularly bisects the vector S −S connecting 1 2 the two foci. It is a remarkable but easily shown property that the projection of this ellipse onto the xy plane is a circle. It follows that the joints are arcs of ellipses and that their projections onto the xy plane constitute a planar diagram, to be called G, whose edges are circular arcs. The diagram G di- vides the plane into cells j each of which is associated with a specific seed S . j A snapshot of such a diagram is shown in Fig.2. All vertices are trivalent, but not all edges end in vertices: some form full circles. Cells may not be convex; they may not be simply connected and may even be disconnected. The projection s = (x ,y ) of S may or may not be in cell j. j j j j Diagrams whose edges are circular arcs (and their higher dimensional generalizations) have been called Möbius diagrams by Boissonnat et al.[13], since they constitute an ensemble that is invariant under Möbius transforma- tions. In the present case where the seeds are Poisson distributed, it seems appropriate to call G a Poisson-Möbius diagram. This diagram is a random object and, since the seed density ρ may be scaled away, it does not depend on any parameter. There are many interesting questions that one may ask about it. 6 3.2 Connection to weighted Voronoi diagrams For a given point r = (x,y) in the xy plane one may ask to which cell j it belongs. This is obviously the cell of seed S whose paraboloid is lower jmin than all the others at r, that is, j = argmin P (r). (3.2) min j j We may rephrase this as a two-dimensional problem in the following way. We refer to the projection s = (x ,y ) of seed S as a two-dimensional seed. i i i i Then r is in the cell of the seed s to which it is closest according to the jmin modified distance function given by1 dist(r,s ) = z−1|r−s |2 +z , (3.3) i i i i inwhichthez arenowinterpreteda‘weights’thatrenderthetwo-dimensional i seed s inequivalent. The diagram G hence appears as an ordinary two- i dimensional Voronoi diagram but with the Euclidean distance replaced by the modified expression (3.3) that weights the seeds. Weighted Voronoidiagramswithavariety ofdistancefunctions have been considered sincemanydecades, oftenmotivatedbypractical applications(see e.g. Ref.[14]). Okabeet al.[1] discuss thestateoftheartofweighted Voronoi diagrams up to the year 2000. Shortly after that, Boissonnat and Karavelas [15] introduced the distance function dist(r,s ) = λ |r −s |2 −µ , (3.4) i i i i i where λ and µ are weights. With this distance definition (3.4) it is easily i i shown that the edge separating the Voronoi cells of two seeds at s and s i j is an arc of a circle. The distance function of this paper, Eq.(3.3), is the special case of Eq.(3.4) with λ = z−1 and µ = −z . i i i i The literature that deals with weighted Voronoi diagrams is often con- cerned either with fairly abstract mathematical properties; or with questions about the computational complexity of algorithms that construct a diagram from a given set of N seeds and their weights. Here we address the subject from a statistical point of view, the diagram G defined above being stochas- tic. Various of its properties may be calculated. Below we will focus directly on the particular property that we need, viz. the total edge length λ per unit area. 3.3 Edge length per unit area in G The question of interest to us here is: what is the total length λ of the edges of G per unit area? For the two-dimensional Poisson-Voronoi diagram of 1 1 A factor in Eq.(3.1) may be ignored without changing the cell structure. 2 7 1/2 seed density ρ the value λ = 2ρ is part of a long list of exactly established 2 2 results. This quantity has the dimension of a length per area, that is, of an inverse length. Hence in the present case we should have λ = cρ1/3 and the nontrivial part of the problem is to calculate the dimensionless coefficient c. Let us consider the infinitesimal line segment connecting (0,0) to (∆x,0). We ask for the probability, to be called p(θ)∆x∆θ, that this line segment be intersected by an edge of G that has an orientation (by which we will mean ananglewith they axis) in [θ,θ+∆θ] (see Fig.3). If we take theline segment to be the side of a parallelogram of height ∆w in the xy plane, we see that the edge length of G inside this parallelogram equals ∆ℓ(θ) = ∆w/cosθ. Therefore the expected edge length crossing a surface area ∆A = ∆x∆w at an angle θ is ∆w h∆ℓ(θ)i = p(θ)∆θ∆x· cosθ p(θ) = ∆θ∆A. (3.5) cosθ The total expected edge length λ∆A crossing a surface element ∆A is equal to the integral of (3.5) on θ, whence π/2 p(θ) λ = dθ . (3.6) cosθ Z−π/2 We will now calculate p(θ). The probability for two arbitrary paraboloids P and P to contribute j k to the surface Γ a joint whose projection crosses the above infinitesimal line segment is also the probability that the joint between P and P crosses the j k strip of zero thickness and infinitesimal width defined by 0 < x < ∆x, y = 0, and z > 0, which in turn is N times the probability that the joint between 2 the paraboloids P and P does so. Let the joint intersect the xz plane in 1 2 (cid:0) (cid:1) (x ,0,z ) and let its projection onto the xy plane intersect the x axis at an 12 12 orientation θ . Hence, averaging over all seed configurations, we have 12 N(N −1) L L p(θ)∆θ∆x = dx dx dy dy dz dz χ (x ,θ ) 1 2 1 2 1 2 12 12 12 2VN Z−L Z0 L N L N N × dx dy dz Θ (z ), (3.7) i i i i 12 Z−L i=3 Z0 i=3 i=3 Y Y Y in which z is the common value of P and P at the point of intersection, 12 1 2 that is z = P (x ,0) for j = 1,2; the Heaviside function Θ imposes that 12 j 12 i 8 θ ∆w ∆l 0 x ∆x 12 y x Figure 3: An element (heavy line) of an edge of the diagram G of circular arcs intersecting a line segment of length ∆x on the x axis. It contributes a length ∆ℓ = ∆w/cosθ to the parallelogrammatic area ∆A= ∆x∆w. the ith paraboloid does not intersect the strip in a point lower than z , that 12 is 1 if P (x ,0) > z , i jk jk Θ (z ) = (3.8) i 12 ( 0 otherwise; and 1 if 0 < x < ∆x and θ < θ < θ+∆θ, 12 12 χ (x ,θ ) = (3.9) 12 12 12 ( 0 otherwise. The N −2 triple integrals on (x ,y ,z ) for i = 3,...,N may be carried out i i i independently for each i. The factor Θ imposes that the integrand vanishes i if (x ,y ,z ) is inside the sphere of radius z around (x ,0,z ). Hence the i i i 12 12 12 result of these integrations is (V − 4πz3 )N−2 which in the limit L → ∞ 3 12 becomes VN−2exp(−4πz3 ρ). Upon taking the limit L → ∞ in (3.7) we get 3 12 ρ2 ∞ ∞ p(θ)∆θ∆x = dx dx dy dy dz dz χ (x ,θ )exp(−4πz3 ρ). 2 1 2 1 2 1 2 12 12 12 3 12 Z−∞ Z0 (3.10) Near the origin of the xy plane the paraboloids P (x,y) with j = 1,2 may j be linearized according to P (x,y) = r +s x+t y +O(x2,y2) (3.11) j j j j where z x2 +y2 x y j j j j j r = 1+ , s = − , t = − . (3.12) j 2 z2 j z j z j j j (cid:16) (cid:17) 9 We now transform from (x ,y ,z ) to new variables of integration (r ,s ,t ), j j j j j j where j = 1,2. The Jacobian is ∂(r ,s ,t ) (1+s2 +t2)3 j j j j j = . (3.13) ∂(x ,y ,z ) 8r2 j j j j With this transformation Eq.(3.10) becomes ∞ ds ds dt dt p(θ)∆θ∆x = 32ρ2 1 2 1 2 (1+s2 +t2)3(1+s2 +t2)3 Z−∞ 1 1 2 2 ∞ × dr dr r2r2χ (x ,θ )exp(−4πz3 ρ) (3.14) 1 2 1 2 12 12 12 3 12 Z0 in which χ couples the variables with indices 1 and 2. The coordinate x 12 12 of the point of intersection is the solution of P (x ,0) = P (x ,0) which 1 12 2 12 upon linearization gives r −r 2 1 x = − . (3.15) 12 s −s 2 1 Using (3.15) we can rewrite the condition 0 < x < ∆x as 12 r −∆x(s −s ) < r < r , s > s , 1 2 1 2 1 2 1 r < r < r +∆x(s −s ), s > s . (3.16) 1 2 1 1 2 1 2 In both cases r is integrated on an infinitesimal interval of length |s − 2 1 s |∆x located at r . This takes account of the condition on x implied by 2 1 12 χ (∆x,∆θ) and shows furthermore that z = r + O(∆x) = r + O(∆x). 12 12 1 2 Hence Eq.(3.14) becomes, after we divide it by ∆x and scale ρ out of the integrand, ∞ |s −s | 1 1 2 ang p(θ)∆θ = 32I ρ3 ds ds dt dt χ (θ ) 0 1 2 1 2 (1+s2 +t2)3(1+s2 +t2)3 12 12 Z−∞ 1 1 2 2 (3.17) in which ∞ I = drr4 exp(−4πr3) = 2(4π)−53Γ(2) (3.18) 0 3 9 3 3 Z0 ang and where χ (θ ) imposes the remaining condition θ < θ < θ+∆θ. The 12 12 12 point of intersection (x ,0,z ) being known, we now look for the line of 12 12 intersection by setting x = x +δx and y = δy. Substituting in (3.11) and 12 eliminating r = r we find that s δx+t δy = s δx+t δy, whence 1 2 1 1 2 2 δx t −t 2 1 tanθ = = − . (3.19) 12 δy s −s 2 1 10