Many physical systems can be studied as collections of particles embedded in space, often evolving in time. Natural questions arise concerning how to characterize these arrangements—are they ordered or disordered? If they are ordered, how are they ordered and what kinds of defects do they possess? Voronoi tessellations, originally introduced to study problems in pure mathematics, have become a powerful and versatile tool for analyzing countless problems in pure and applied physics. We explain the basics of Voronoi tessellations and the shapes that they produce and describe how they can be used to characterize many physical systems.
I. INTRODUCTION
Physical systems can often be abstracted as large sets of pointlike objects in space, with each point representing, for example, the position of an atom, colloidal particle, organism, or celestial body. Fundamental questions arise when studying these systems regarding the manner in which their constituent particles are organized. Are they ordered or disordered? If they are ordered, are the particles arranged in a crystallike fashion, and if so, what kind? Systems that are primarily ordered often contain defects of varying dimension and kind, all of which can impact the macroscopic static and dynamic properties of a system. How can these defects be classified? In systems that are nominally disordered, such as liquids, glasses, and granular materials, how can the disorder be described in a practical and effective manner?
Voronoi tessellations provide a natural tool for converting questions about the arrangements of particles into ones about polygons and polyhedra, possibly irregular, and about how they fit together to tile space. The analysis of these shapes can tell us much about particlelike systems on a wide range of scales.
II. VORONOI CELLS
In many applications, the space M under consideration is $ \mathbb{R} 2 $ or $ \mathbb{R} 3 $, or some subregion thereof, and the metric d is the standard Euclidean one. In two dimensions, Voronoi cells are convex polygons that together partition space, while in three dimensions, they are convex polyhedra.
Figure 1 shows examples of particles and their Voronoi cells in two and threedimensional systems. The Voronoi cells have different areas and volumes as well as different numbers of edges and faces. The shape of each Voronoi cell describes the order in the vicinity of a particle and in the system more generally. In perfect crystals, for example, all Voronoi cells are identical, whereas in disordered ones, there can be significant variation of their properties.
Voronoi cells also provide a scalefree mechanism to characterize adjacencies between particles: two particles are defined to be neighbors if and only if they share a Voronoi edge (in $ \mathbb{R} 2 $) or face (in $ \mathbb{R} 3 $). If edges are drawn between all such neighbors, this results in the Delaunay triangulation, another widely used construction in computational geometry^{2} (see Sec. IV E). Voronoi tessellations and Delaunay triangulations are sometimes referred to as dual objects, because each can be used to construct the other. Decomposing a region into polyhedral tiles through the Voronoi tessellation has proven to be a powerful tool in many branches of science and engineering.^{3}
A. Construction
The definition given by Eq. (1) does not directly lend itself to a practical algorithm for constructing Voronoi cells. It is, thus, productive to also think about Voronoi cells from an alternative perspective. Consider the two particles in Fig. 2(b) and the line separating them. On one side of the line lies the region of space closer to the magenta particle, and on the other side lies the region of space closer to the yellow one; points on the line itself are equidistant to the two particles. If these were the only particles in the system, we would be left with two Voronoi cells, one on each side of the dividing line. Each of these Voronoi cells is called a halfspace. In threedimensional spaces, the two Voronoi cells are separated by a plane instead of a line.
Particles located far away from a central particle do not influence its Voronoi cell, because associated halfspaces completely contain the intersection of halfspaces associated with closer particles. The shape of a Voronoi cell is, thus, generally determined only by nearby particles.
The definition provided by Eqs. (4) and (5), thus, motivates a simple constructive algorithm for building Voronoi cells: iteratively cut the region around a particle s_{i} by intersecting halfspaces. Figure 3 illustrates the cutting of a threedimensional Voronoi cell through a halfspace associated with an added particle.
B. Algorithmic and computational aspects
The time and memory necessary to compute and represent a Voronoi tessellation of an entire system increases with the number of particles. In a system of n particles, constructing the Voronoi cell of each particle requires, in theory, computing the intersection of n − 1 halfspaces, and so calculating all Voronoi cells would require consideration of $ O ( n 2 ) $ halfspaces.^{10} However, as described in Sec. II A, usually only the halfspaces for nearby particles influence the Voronoi cell. This consideration allows efficient algorithms to be designed^{11} that first initialize the candidate Voronoi cell to fill the entire domain, and then sweep outward, performing halfspace intersections for particles of increasing separation, until a termination criterion is reached and it is possible to guarantee that any particles further away cannot possibly affect the Voronoi cell shape. For example, if the candidate Voronoi cell is contained within a sphere S of radius R, and all particles within the sphere of radius 2R have been considered, then we can guarantee that the Voronoi cell computation is complete, because the halfspace for any particle further away will completely contain S.
For typical particle arrangements that are evenly distributed, each Voronoi cell requires O(1) halfspace intersections, meaning that the total computation time scales linearly with the number of particles. However, in certain pathological cases, the computation time may not be linear. For example, in three dimensions, define $ L = { \u2212 n , \u2212 ( n \u2212 1 ) , \u2026 , n} $, and consider the two sets of particles $ P 1 = { ( 1 , k , 0 ) \u2009  \u2009 k \u2208 L} $ and $ P 2 = { ( \u2212 1 , 0 , k ) \u2009  \u2009 k \u2208 L} $. Then, the Voronoi cell for a particle in P_{1} contains a face for every particle in P_{2} and vice versa, and, thus, each Voronoi cell requires O(n) halfspace intersections in this case, meaning that the total computation time will scale at least quadratically.
In two dimensions, Voronoi cells can also be constructed using a “walking method.”^{12,13} This approach traces around the edges and vertices of the Voronoi cell, directly constructing it without explicitly considering the halfspaces. Efficient implementations allow this procedure to be performed in O(1) time for a typical Voronoi cell.^{14} However, unlike halfspace intersections, this method does not have an obvious generalization to higher dimensions.
Because each Voronoi cell can be computed and analyzed independently, the cellbased perspective is wellsuited to current computational hardware, where parallelism is emphasized. A typical workflow could involve calculating a Voronoi cell, storing its statistics, and then deleting it and moving onto the next. On a multicore or parallel computer, each thread or process can compute subsets of Voronoi cells independently of the rest. A drawback of the cellbased approach is that it lacks explicit topological information about how the cells are connected together—it is not straightforward to identify common faces, edges, and vertices between cells. This can be a problem once the effects of numerical rounding error of floating point arithmetic^{15} are accounted for. In a typical situation in two dimensions, three Voronoi cells will meet at a single vertex [see Fig. 4(a)]. However, if the particles are precisely aligned, then it is possible that four or more cells meet at one vertex. In this case, due to numerical rounding error, one or more Voronoi cells may have additional small faces, leading to inconsistent topology between the cells [see Fig. 4(b)]. This problem is not unique to $ \mathbb{R} 2 $ and happens in arbitrary dimensions. In $ \mathbb{R} n $, it occurs when n + 1 precisely aligned particles are equidistant to a single point.
There are computational methods for eliminating these discrepancies. For example, the Voronoi cells can be computed using an exact arithmetic system (e.g., GMP^{16}) that eliminates rounding error, albeit for substantially higher computational cost. However, in many cases, this problem has a limited effect on the results of interest. Cellbased computations, such as surface area and volume, are not sensitive to small topological changes. Furthermore, even though these numerical errors may introduce spurious faces and edges, it is still possible to develop a topological framework in which such errors are anticipated and accounted for (see Secs. III B and III C).
Alternatively, algorithms exist that can compute the full Voronoi tessellation as a single mesh, i.e., returning a full graph of blue edges as shown in Fig. 1. These approaches require building a global data structure for the mesh and, thus, are less amenable to parallelization when compared to the cellbased perspective. When the Voronoi cells are extracted from this mesh they will have internally consistent topology. Some popular methods for computing the full Voronoi tessellation include the Fortune sweeping algorithm^{17,18} and the incremental approach, whereby the mesh is continually updated as new particles are added.^{19,20} Another method is to use liftup mapping, projecting a particle $ s \u2208 \mathbb{R} n $ to a paraboloidal surface $ ( s ,   s   2 ) \u2208 \mathbb{R} n + 1 $. The hyperplanes tangential to the surface form facets exactly matching the Voronoi tessellation when projected back to $ \mathbb{R} n $. This method can be efficiently implemented in arbitrary dimensions using the Quickhull algorithm.^{21}
A variety of software packages are available for computing the Voronoi tessellation. The Qhull library^{22} computes Voronoi cells using the Quickhull algorithm and forms the basis of Voronoi computations in MATLAB. The Computational Geometry Algorithms Library (CGAL)^{23} has routines for calculating a variety of geometrical constructions, including the Voronoi tessellation. Voro++ is a software library for performing cellbased computations and has both a commandline utility and a C++ application programming interface for calculating statistics of Voronoi cells.^{24,25} Voro++ is incorporated into the OVITO software package for visualizing and analyzing particle systems^{26} and is integrated into the Large Atomic/Molecular Massively Parallel Simulator (LAMMPS) for particle simulation.^{27–29} Several parallel implementations are also available. Starinshak et al.^{30} use a parallel divideandconquer approach, where each processor computes Voronoi cells in a subdomain, after which they are joined together. Other examples include PARAVT^{31} and a GPU implementation.^{32}
III. THE SHAPES OF VORONOI CELLS
As the Voronoi cell of a particle is constructed using information about neighboring particles, its geometric and topological features, in turn, reflect important structural information about the spatial arrangements of particles, both locally and globally. We consider several examples.
A. Geometric aspects
In continuum mechanics, density is modeled as a property that varies continuously in space. This viewpoint is difficult to justify when modeling systems at the atomic scale, because mass in such systems is localized at discrete points. A reasonable question, thus, arises as to how density should be defined and analyzed, and the volume of a Voronoi cell is often used for this purpose. In finite systems, the average area or volume of all Voronoi cells is equal to one divided by the average particle density. The variation of particle density throughout a system can also be captured through the geometric properties of Voronoi cells. Voronoi cells with small areas and volumes indicate particles in densely packed regions; conversely, those with large areas or volumes indicate particles whose neighbors are located far away. Section IV C discusses an example of using Voronoi volumes to analyze a granular material simulation.
The variation in the areas and volumes of Voronoi cells can also be indicative of structural abnormalities. In lowtemperature crystalline systems, for example, the Voronoi cells of bulk crystal particles are nearly identical. Therefore, Voronoi cells with significantly larger or smaller areas and volumes can be indicative of localized defects such as vacancies and interstitials as well as dislocations and grain boundaries.^{33}
In addition to the geometric features of individual Voronoi cells, the consideration of their distributions can also reflect information about global order. In perfect crystals, for example, all Voronoi cells have identical volumes and surface areas. At low temperatures, these distributions are narrowly concentrated around their mean values. As the system is heated and becomes increasingly disordered, these distributions become broader. These differences can be observed in Fig. 5, which illustrates both ordered and disordered systems. Although the Voronoi cell areas of particles in the crystals are nearly identical, there is considerable variance in those belonging to disordered regions.
B. Topological aspects
Topology is the mathematical study of connectivity and shape,^{34} and a thorough analysis of the topologies of individual Voronoi cells, as well as their distributions over entire systems, can provide much insight about local arrangements and global ordering. Unlike geometric properties such as area and volume, topological ones, such as the number of edges or faces of a Voronoi cell, are typically discrete, even nonnumeric. Indeed, the manner in which faces of a threedimensional Voronoi cell are arranged relative to one another is captured through the isomorphism class of the cell's edge graph; this classification is typically represented in the language of graphs instead of numbers.^{35}
In solidstate physics, the unique Voronoi cell of a crystal lattice is often referred to as the Wigner–Seitz cell of the system.^{36,37} In some twodimensional crystals, for example, every Voronoi cell is a hexagon. Particles whose Voronoi cells have more than or fewer than six edges can be identified as belonging to defects. For example, Figs. 5(a)–5(c) illustrate crystalline systems in which most Voronoi cells have six edges. In contrast, Voronoi cells of many particles lying along the grain boundary in Fig. 5(b) have more or fewer edges.
In contrast to crystals, nominally disordered systems such as liquids and gases exhibit a substantially wilder set of particle arrangements, reflected in a larger variation of their topological features. Instead of a single polyhedralshaped tile centered at each lattice point (i.e., the Wigner–Seitz cell), disordered systems are tiled by a large assortment of Voronoi polyhedra. All simple polyhedra, for example, appear with some frequency as Voronoi cells in the ideal gas,^{38} providing an endless pool of Voronoi cell types for classifying particle arrangements.
The topology of a Voronoi cell in two dimensions can be characterized by its number of edges. This classification can be refined by further considering the number of edges of the Voronoi cell's neighbors; such a classification is adopted in Fig. 5. The topology of a Voronoi cell in three dimensions is more delicate due to the added complexity of describing the manner in which faces of a polyhedra are arranged. This task can be done using an algorithm of Weinberg^{39} originally designed to determine whether two planar graphs are isomorphic. Because edge graphs of threedimensional Voronoi cells are always planar,^{40} Weinberg's graphtracing algorithm provides a systematic description of Voronoi cell topologies.^{35,41} This classification of local structure avoids some inherent limitations of more conventional methods.^{41,42} Furthermore, it facilitates the a priori enumeration of Voronoi topologies associated with particle structures.^{41}
The VoroTop^{43} program is a commandline utility for calculating Voronoi topologies of threedimensional particle systems. The program uses the Voro++ software library^{24} to compute the individual Voronoi cells and then performs additional topological analysis. By considering special families of Voronoi topologies associated with given crystalline systems, VoroTop facilitates the automated identification of defects in crystalline systems as well as automated statistical characterization of disordered systems.
C. Stability
If the properties of Voronoi cells are to be useful for classifying local structure, then it is important to understand the manner in which those properties can change under small perturbations of particle positions. This issue is important to consider because measurement error is inevitable when working with experimental data. Even numerical codes, in which particle coordinates are known precisely, are subject to rounding and approximation error.
The volume of a Voronoi cell is a continuous function of particle positions, and small perturbations of the particle coordinates, thus, result in small changes in cell volumes. Likewise, most geometric features of Voronoi cells, such as perimeters and surface areas, are continuous functions of particle positions.^{44} When geometric quantities are used for classification purposes, thresholds for the range of possible values must be chosen. Even small errors in the measurements of particle positions can result in classification errors.
Unlike geometric properties of Voronoi cells, topological ones can change abruptly under small perturbations of particle coordinates.^{45,46} For example, the number of edges or faces of a Voronoi cell can change under arbitrarily small perturbations of particle positions. Figure 6 illustrates a perfect and perturbed square lattice in which particles are colored according to the number of edges of their Voronoi cells. Small perturbations of the particles, which can be seen as representative of thermal vibrations or measurement errors, result in Voronoi cells that are topologically distinct from those in the unperturbed case.^{47,48} This instability arises from symmetries of the lattice which result in points equidistant to four neighboring particles. Small perturbations break this symmetry, resulting in topological changes.
This instability can be understood by considering an abstract configuration space of all possible local arrangements of particles,^{41} schematically illustrated in Fig. 7. Voronoi topology provides a natural means of segmenting this space into regions representing particles with a given local structure. Particles in unstable crystals correspond to points located on the boundary of several such regions, and so perturbations of particle positions, corresponding to small perturbations in this space, typically result in different Voronoi cell topologies. The set of all topologies that can result from small perturbations can then be associated with a given structure type, allowing for a simple method for dealing with experimental error and thermal vibrations.^{41} Combinatorial factors limit the number of types that can result from small perturbations, corresponding to the number of regions meeting at the representative point in the configuration space. In many cases, a complete set of types that can result from small perturbations can be determined analytically or numerically.^{41,49}
Not all lattices and crystals are unstable in this sense. In two dimensions, for example, the hexagonal lattice, whose unique Voronoi cell is a regular hexagon, is stable under small perturbations; cf. Fig. 5(a). This arrangement of particles corresponds to a point in the interior of a region of configuration space with fixed Voronoi topology. In three dimensions, the bodycentered cubic lattice is stable, whereas the majority of common crystals, such as the simple cubic, facecentered cubic, hexagonal close packed, and diamond crystals, are not.^{41}
IV. APPLICATIONS
In recent decades, Voronoi tessellations and their generalizations have been used to analyze countless problems throughout the physical sciences, as well as in the life sciences, engineering, and pure mathematics. For example, Voronoi tessellations have found broad use in computational fluid dynamics,^{50–52} secure data compression and transmission,^{53–56} machine learning,^{57–59} quantum error correction,^{60,61} and a wide range of biological, chemical, and ecological problems.^{62–65} Even a brief survey of the dozens of such applications is well beyond the scope of this article. Instead, we mention several examples of particular interest to physicists and detail how Voronoi tessellations are used in their analysis.
A. Crystals and defects
Many common materials such as diamonds, table salt, and most everyday metals, are crystalline in nature. Although the bulk of these materials are highly ordered, the underlying structure is punctuated by defects of various dimensions and kinds.^{37,66} Material scientists wish to understand these defects, because they impact macroscopic material properties such as strength and ductility.^{67} To study defects, computational material scientists often use supercomputers to simulate systems with millions, or even billions, of atoms. However, without a way to automate the accurate identification of the defects, the massive amount of data is not useful. Automating the analysis of large atomistic data sets, thus, motivates the development of methods to identify and classify defects using local structural data.
Figure 8 illustrates a twodimensional system with two crystals, a boundary between them, and a localized vacancy. Particles are colored according to their Voronoi cell areas and topologies in Figs. 8(b) and 8(c), respectively. Although the vacancy defect is associated with Voronoi cells with aboveaverage areas, corresponding to a decrease in the local density, area changes resulting from thermal fluctuations of other particles in the system make it difficult to identify the defects by the analysis of Voronoi cell areas alone. When colored according to Voronoi topologies, particles belonging to the two bulk crystals appear locally ordered, and their Voronoi cells all have six edges. Those adjacent to the grain boundary and vacancy have more than or fewer than six edges. This observation suggests using the topology of Voronoi cells to identify, and possibly classify, defects in crystalline systems.
Although this example is twodimensional, this approach is also effective in three dimensions and enables the accurate identification of defects in crystalline systems at high temperatures, where conventional methods fail.^{41,43,49}
B. Disordered systems
Although originally introduced to study crystalline systems, Voronoi cells also have a long and storied history in the study of liquids, glasses, and other nominally disordered systems. In the early 1930s, John D. Bernal, the great pioneer of xray diffraction and founder of protein crystallography, initiated a study of simple liquids.^{68–70} Bernal suggested looking at the Voronoi cells of particles, and considering not only their number of faces but also the number of edges of each face. This detailed analysis of Voronoi cells enabled him to describe the structure of liquids through a topological description of their constituent particles.
Despite lacking crystalline order, nominally disordered systems can still be distinguished from one another and analyzed through local structural features. For example, liquid copper and an ideal gas are both disordered, and yet are structurally distinct. One of the challenges in studying disordered systems is meaningfully describing and accurately identifying different kinds of structural order.
The distribution of Voronoi topologies can be effectively used for this purpose.^{71–74} In particular, a detailed list of the Voronoi topologies and their observed frequencies provides important structural insight into local structure character. Icosahedral symmetry, studied in the context of liquids and glasses,^{75} is nothing but a single type of local structure, and has shed light on the tension between local energyminimizing and entropymaximizing principles in determining how particles arrange in disordered systems. Although our understanding of the distribution of local Voronoi topologies in general systems is still incomplete, these data have been recently used to identify phase transitions in supercritical fluids.^{76–78}
C. Granular flows
Voronoi tessellations have also been successful in characterizing the dynamics of granular flows. Despite being common in everyday experience, granular materials have surprisingly complex mechanical properties: in some situations, they behave like solids and support stress, whereas in other situations, they can flow like liquids.^{79} This complex behavior has prompted much study in the past several decades to fully understand the mechanics of granular media.^{80,81} A key quantity of interest is how closely packed the particles are. If they are tightly packed, they will become jammed^{82} and behave like a solid, but in order to flow, there must be some free volume for the particles to move past each other. The notion of free volume has been extensively studied for many amorphous materials^{83,84} and has also been used as the basis for simple models of granular drainage.^{85,86}
Voronoi cells are useful for precisely characterizing the local volume fraction, i.e., what proportion of space is occupied by the particles. A simple approach is shown in Fig. 9(a), where a small pink square of volume V_{s} is drawn. Twelve particles with volume V_{p} have centers inside the square. Thus, the local volume fraction can be estimated as $ 12 V p / V s $. However, this measurement is sensitive to small changes in the particle positions: if the center of a particle moves outside the box, the volume fraction will become $ 11 V p / V s $, which is a 8% relative change. This change is problematic, because small fluctuations of the packing fraction of several percent can have a large effect on the mechanical properties. This calculation can be improved by computing the Voronoi cells for the particles. The pink square volume can be replaced by the total volume V_{v} of the Voronoi cells [Fig. 9(b)], resulting in a more faithful representation of the space associated with the particles and a more accurate computation of the local volume fraction.
To illustrate this method, we performed a granular discreteelement simulation, where the position, velocity, and angular velocity of each particle are integrated using Newton's laws. We simulate 140 000 particles of diameter d using LAMMPS,^{27,28} where parameters in the interparticle contact model are derived from Rycroft et al.^{87} with a Coulomb friction coefficient of 0.5. Particles are first poured into a container of width 140d and thickness 12d using periodic boundary conditions in the thickness direction. Figure 9(c) shows a snapshot of the simulation, where the bottom layer of particles shown in yellow are frozen in place. The remaining particles are all physically identical but are colored in alternating bands of width 8d. Figure 9(d) shows the local packing fraction, achieving a value of $ \u2248 63 % $, which is typical of random close packing.^{88}
A section of the frozen particles are then slowly raised a distance of $ 4.7 d $. Figure 9(e) shows how the layers of particles are deformed, but it is visually difficult to see any change in the packing fraction. However, the Voronoibased calculation demonstrates that two bands of lower packing fraction form in regions undergoing high shear, where particles must move apart slightly in order to move past one another. Such analyses can be used to probe the mechanics of granular media.^{89}
Voronoi volumes have been used to probe a variety of aspects of granular media. Although Fig. 9(b) demonstrates the computation of the local packing fraction in a small region, the same principle can be used to define the local packing fraction at the level of each particle by dividing the particle volume by its Voronoi cell volume. This approach has been used to examine local volume fluctuations^{90–93} and dynamic jamming fronts.^{94} The use of Voronoi cell faces to identify neighboring particles^{95} and the examination of Voronoi cell anisotropy^{96,97} have also proven fruitful in the study of granular media and other amorphous systems.
D. Astronomy
In the examples we have considered so far, the particles around which the Voronoi cells are constructed have represented atoms or other smallscale objects. Voronoi tessellations have also found applications when the particles under consideration are of astronomical scale. In an attempt to understand the Universe and its origins, astronomers have spent recent decades imaging and cataloging the observable sky, carefully mapping out its stars, galaxies, and other celestial bodies. Dozens of digital surveys, such as the Two Micron All Sky Survey^{98} and the Sloan Digital Sky Survey,^{99} have generated many petabytes and exabytes of data from which scientists have culled novel insight into many secrets of the cosmos.
To fully exploit the breathtakingly large amount of recorded data, scientists have been eager to develop automated computational techniques to aid in its analysis. Voronoi tessellations have consequently been used to study the large scale distribution of galaxies^{100} as well as to identify individual astronomical objects such as galaxy groups and clusters,^{101} haloes,^{102} voids,^{103} and filaments.^{104} Each of these studies has provided insight into the early history of the universe, its evolution, and its largescale structure.
E. Numerical analysis
One important problem that arises in modeling countless physical systems is the numerical solution of partial differential equations on domains with nontrivial shapes. A necessary first step in solving many such problems is the discretization of the domain, and Voronoi tessellations are commonly used for this purpose.^{105–108} Here, we illustrate how to solve one such problem using both a uniform mesh and an adaptive one, both constructed using Voronoi tessellations.
We aim to solve Laplace's equation where the Dirichlet condition is set to $ u ref $ on the boundary $ \u2202 \Omega $. Because the solution to the Laplace equation is unique in this case, it must equal $ u ref $ everywhere on Ω. Thus, we can compare our numerical solution $ u \u0302 ( x , y ) $ to $ u ref $ to test the accuracy of the numerical method. We use a piecewise linear triangular mesh.^{109} To generate such a mesh, we first create a Voronoi tessellation and then use it to construct the dual Delaunay triangulation by connecting points that share a Voronoi edge.^{2} The Delaunay triangulation has the property that for any triangle, no other points are in the circumcircle of the triangle. Delaunay triangulations are attractive for use with finite element methods, because the definition disfavors thin triangles that have small internal angles and lead to large numerical errors.
We, therefore, need to compute a Voronoi tessellation of Ω using evenly distributed points. To do this, we use a centroidal Voronoi diagram,^{112} where each point coincides with the weighted center of mass, or centroid, of its Voronoi cell. A centroidal Voronoi diagram can be computed via Lloyd's algorithm.^{112,113} For an initial arbitrary set of points, Lloyd's algorithm iteratively moves the points to transform the corresponding Voronoi diagram into a centroidal Voronoi diagram. The algorithm is as follows:

Choose n initial points and compute their associated Voronoi diagram.

Compute the centroid for each of the Voronoi cells.

Move each point to its Voronoi cell's weighted centroid location.

Repeat steps (2) and (3) until a stopping criterion is satisfied, such as a maximum number of iterations or a minimum distance of point movement.
Figure 11 illustrates how Lloyd's algorithm works. As seen in Fig. 11(a), we first choose 11 random points in the unit square and compute their corresponding Voronoi diagram. We then implement Lloyd's algorithm for 20 iterations. After 20 iterations, as seen in Fig. 11(d), we obtain a Voronoi diagram where the points are uniformly distributed in the domain, and their Voronoi cells are of similar sizes.
Figure 10(a) shows a centroidal Voronoi diagram of Ω using n = 1468 points. From here, the corresponding Delaunay triangulation can be computed, and the finite element method can be implemented. The numerical solution $ u \u0302 $ is shown in Fig. 10(b). Figure 10(c) shows the absolute errors, $  u \u0302 ( x , y ) \u2212 u ref ( x , y )  $, for the numerical solution, indicating that the largest errors are near the four concave creases. This is expected because the solution varies most rapidly there.
The density field can be set manually, but here we use an automatic procedure that forces the triangles in the interior of the domain to be large, the triangles near the boundary to be small, while ensuring a smooth gradation in mesh size.^{106,114–116} Figure 10(d) shows the resulting Voronoi tessellation using the same number of n = 1468 points. The points of the resulting CVD distribute more densely around the corners and the concave creases of the shape and less densely in the interior of the three lobes. Figure 10(e) shows the corresponding Delaunay triangulation and the finiteelement solution. Because the mesh resolution is increased in the regions where $ u ref $ varies most rapidly, the numerical errors are reduced when compared to the uniform mesh, as shown in Fig. 10(f).
V. CONCLUSIONS
Voronoi tessellations provide a natural framework for describing and analyzing arrangements of particles on a wide range of scales. Careful analysis of geometric and topological features of their constituent polyhedronshaped cells can provide insight not readily apparent from other perspectives.
Despite being introduced over a century ago as an abstract mathematical concept, Voronoi tessellation continues to find new applications in many contemporary research problems. In the last decade, there has been much interest in datadriven approaches in science, such as the analysis of large databases,^{117} and the use of new machine learning methods.^{118} These approaches often rely on the construction of simple descriptors of a system, and Voronoi cell features have proven highly useful, such as for screening highthroughput databases^{119,120} and for predicting properties of crystalline structures.^{121}
VI. SUGGESTED PROBLEMS
We suggest several problems suitable for motivated students, perhaps as part of a course in computational physics. Those marked with a $ \u22c6 $ might form the basis of a semesterlong research project.

Accurately measuring particle positions in experimental systems can be difficult. Drawing quantitative conclusions based on Voronoi tessellations of such data, thus, presents significant challenges. Let's anticipate the error of Voronoi cell features as a function of measurement error in a simple system. Consider a square lattice in the plane, and displace each of its points by adding an independent and identically distributed random perturbation, such as a twodimensional Gaussian. How do the distributions of areas and perimeters depend on the magnitude of the perturbations?

Repeat problem (1) in three dimensions for a cubic lattice and consider the distributions of cell volumes and surface areas. How do they depend on the magnitude of the random perturbations? How do the distributions of topological features of the Voronoi cells, such as the number of edges or faces, depend on the magnitude of the perturbations?

Images of the Voronoi cells can be computed with a simple procedure in two dimensions. Consider a number of particles $ { s 1 , s 2 , \u2026} \u2208 \mathbb{R} 2 $, and define an image covering a rectangle in $ \mathbb{R} 2 $. For each pixel location x in the image, assign it a different color based on computing $ arg min j d ( x , s j ) $ to find the nearest particle.^{122} This procedure will work for any distance metric. Try computing the Voronoi cells for different pnorms,^{15} where $ d ( x , y ) =   x \u2212 y   p $, and compare their shapes to Voronoi cells with the Euclidean norm.^{123}

Instead of coloring pixels according to the nearest particle as in problem (3), try coloring particles according to the second nearest particle. What do the images look like in this case?

$ \u22c6 $ Computing Voronoi cells in arbitrary domains has several subtleties. Consider the particle arrangement in Fig. 12(a). If the Voronoi cells are restricted to the domain U shown, then one of the Voronoi cells breaks into two disconnected components [see Fig. 12(b)]. Can criteria be placed on the shape of the domain and/or the positions of the particles to avoid this? Can the Voronoi cell definition be modified to ensure that each Voronoi cell is a single connected component?

$ \u22c6 $ Use a molecular dynamics simulator package, such as LAMMPS, to model a crystalline system at temperatures just below melting. What features of the Voronoi cells indicate that the system is still crystalline? Then, heat the system until it has melted. What features of the Voronoi cells reflect the fact that the system has melted?

$ \u22c6 $ Although icosahedral symmetry has been primarily studied in the context of many glasses, symmetric particle arrangements also appear in many other nominally disordered systems. Even an ideal gas, for example, contains Voronoi cells that are very symmetric, including tetrahedra, cubes, and icosahedra. Model a liquid at different temperatures and pressures. What kinds of ordered patterns, indicated by symmetric Voronoi cells, can be observed in these otherwise disordered systems?

$ \u22c6 $ Xray diffraction patterns are often used to characterize ordered and disordered systems. What relations can be found between xray diffraction patterns and the distribution of Voronoi cell features such as areas, volumes, and topologies?
ACKNOWLEDGMENTS
This research was supported by a grant from the United States—Israel Binational Science Foundation (BSF), Jerusalem, Israel through Grant No. 2018/170. C. H. Rycroft was partially supported by the Applied Mathematics Program of the U.S. DOE Office of Science Advanced Scientific Computing Research under Contract No. DEAC0205CH11231.
AUTHOR DECLARATIONS
Conflict of Interest
The authors declare that they have no conflicts of interests related to this article.