Dynamical arrest of topological defects in 2D hyperuniform disk packings

We investigate collective motions of points in 2D systems, orchestrated by Lloyd algorithm. The algorithm iteratively updates a system by minimising the total quantizer energy of the Voronoi landscape of the system. As a result of a tradeoff between energy minimisation and geometric frustration, we find that optimised systems exhibit a defective landscape along the process, where strands of 5and 7-coordinated dislocations are embedded in the hexatic phase. In particular, dipole defects, each of which is the simplest possible pair of a pentagon and a heptagon, come into the picture of dynamical arrest, as the system freezes down to a disordered hyperuniform state. Moreover, we explore the packing fractions of 2D disk packings associated to the obtained hyperuniform systems by considering the maximum inscribed disks in their Voronoi cells.


Introduction
Disordered hyperuniformity describes an amorphous state of matter that suppresses variations in the density of particles on large length scales. Hyperuniform patterns can be readily observed in an extensive range of systems at various scales in nature [1][2][3][4][5][6][7][8]. Recently, hyperuniformity has been utilised as a framework to describe glassy states of jammed granular materials [2-4, 6, 8] giving us a deeper understanding of their out-of-equilibrium structures while being maximally disordered yet mechanically rigid.
In particular, it has been shown that a wide class of maximally disordered jammed packings of hard-particles (in 2D and 3D), including spherical and non-spherical particles (or their 2D counterparts, i.e. disks, ellipses and superdisks) with polydispersity in size, is hyperuniform. That is, the two-point correlation function of these systems decay asymptotically with scaling r −(d+1) , where d is the Euclidean dimension of space, known as 'quasilong-range' correlations [3,4]. Such correlations are contrasted with typical disordered systems in which pair correlations decay exponentially fast. Moreover, in dynamical systems, hyperuniformity has been observed in particle suspensions, where the application of periodic shear can suppress fluctuations in particle concentration across all length scales leading to a hyperuniform state [9].
Hyperuniform point patterns can be constructed using Lloyd's centroidal Voronoi diagram algorithm (Lloyd algorithm, hereafter) in both 2D and 3D with periodic boundary conditions [10]. The algorithm is an iterative process that homogenises a spatial distribution of points as follows: for a given point configuration at step t ≥ 0, the algorithm considers the corresponding Voronoi decomposition and computes the centroid of each Voronoi cell to replace the generating point of the cell. By repeating the process, the given configuration converges to a certain universal disordered state, namely a hyperuniform state. The algorithm optimises the space coverage by maximising the distance between points, while at the same time, minimising quantizer energy associated with the Voronoi cells constructed around each point. This competition results in hyperuniform point patterns, where the system becomes maximally jammed with quasi-long-range correlations.
In this proceedings, we show that 2D hyperuniform configurations obtained by Lloyd algorithm possess defective landscapes, where 5-or 7-coordinated topological defects take part in settling down the competition between energy minimisation and geometric frustration on a flat surface. With the aim of unveiling structural motifs that mediate point-point interactions and drive a dynamical arrest into a defective landscape along the course of quantizer energy minimisation of 2D point patterns, we highlight that dipole defects, each of which is a single pair of a pentagon and a heptagon, exhibit higher number density than other defect components of different sizes. We further investigate the 2D disk packing that arises from these defective hyperuniform structures and compare "optimal hyperuniform" and "ideal hexagonal" disk packings.

Methods
As input to Lloyd algorithm, we generate 20 realisations of Poisson point processes (Poisson-PP) as a statistical en- semble, each with N = 10 4 points in a square box of side length L = 100 (for the unit number density) with periodic boundary conditions. The algorithm runs for at least step 10 4 , where the total quantizer energy of a system remains unchanged with further iterations. The quantizer energy minimisation throughout the algorithm is done by using Papaya2, a computational tool for Minkowski tensors of 2D objects, hence to compute the quantizer energy of Voronoi cells [11]. To investigate the collective reconfiguration of points in a planar random point process into a disordered hyperuniform state, we adopt Crocker-Grier algorithm for a particle tracking velocimetry (PTV) analysis [12]. To obtain the packing fractions of converged hyperuniform systems, we compute the maximum inscribed disk of each Voronoi polygon via a 'shrinking edge algorithm' implemented in MATLAB [13].

Results
The global minimum (quantizer) energy landscape for a 2D system can be achieved by a hexagonal lattice structure with every point having six neighbors around, which is equivalent to the densest packing of hard-spheres in 2D. However, despite the energetic stability of hexagons, we observe along the Lloyd process a dynamical arrest of 5and 7-coordinated topological defects, tending to pair up and create strand-like connected components. These topological defects are crucial to evolve a system into a defective phase by hindering the formation of a perfect hexagonal lattice in the global minimum scheme for quantizer energy and contribute to disordered hyperuniformity. Such a phase is classified under hexatic phase [14]. Each topological defect embedded in a 2D hexatic phase can be characterised by its topological charge 6 − c, where c is the coordination number of the defect point [14]. In view of the Voronoi cell of a defect, we can associate a pentagonal cell with +1 charge and a heptagonal cell with −1. We denote by k the size (number of member defects) of a connected component of pentagonal and heptagonal cells. Hence, if we consider a connected component of three defective cells, i.e. a defect component of k = 3, comprising of two pentagons and one heptagon, its net topological charge is (+1) × 2 + (−1) = +1, while the one with two heptagons and a pentagon has the negative unity charge.
We consider the first step at which the number of pentagonal cells and that of heptagonal cells balance out with each other, so that the net topological charge of a system becomes zero (neutral charge) and denote it by t 0 . For some realisations, the system-wide net charge goes through a small perturbation around zero after step t 0 , but is eventually stabilized to stay zero for the rest of the course. The mean value of step t 0 over the prepared 20 realisations of Poisson-PP is 85.7±18.68.

Dynamical arrest in Lloyd algorithm
We analyse the dynamics of point reconfiguration along the Lloyd process via PTV analysis. In the top panel of Figure 1 we have trajectory plots of the entire system within step ranges (a) [1,100] and (b) [7000, 8000] of the Lloyd iterations. At the early stage of the process, a fast dynamics is prominent system-wide, which becomes slower with iteration time. As the system starts to be jammed into a hyperuniform state at the later stage of the algorithm, much slower dynamics can be observed, leaving only a few localised parts of the system to continue to evolve as in Figure 1(b).
In Figure 1(c), the slowing dynamics is captured by the decaying variance of displacements ∆ t,t+1 over iteration step t, reminiscent of the glassy slow relaxation. The following plot in Figure 1(d), the mean-square displacements (MSD) of points, averaged over 20 different realisations of Poisson-PP, complements the story of the dynamical arrest of points: after passing a rapid dynamical regime at the early stage of evolution, defective cells of pentagons and heptagons start to stabilise the system in a hexatic phase by balancing out the number of one another, resulting in suppressed displacements. Then it is followed by a plateau regime as the iteration step reaches 10 4 , where a system is developed into a hyperuniform configuration. The MSD regime of our systems can be found analogous to relaxation dynamics in glassy systems, fitted by a stretched exponential law, also known as Kohlrausch-Williams-Watts (KWW) law [15][16][17]:

Defective landscapes in hyperuniform systems
Given a Poisson-PP, we show in Figure 2(a) Voronoi landscapes of the evolving point configuration obtained at three selected steps t = 10, t 0 and 10 4 . At an early stage of the process, the Voronoi cells are comprised of a variety of n-gons for n ≥ 3. However, as the points set to arrange into a hyperuniform pattern along the algorithm, Voronoi landscapes are left with n-gons for n = 5 (blue), 6 (white) and 7 (red) only from step t = t 0 , contributing to the neutral charge of the whole system. By the time when the algorithm reaches step 10 4 , pentagonal-and heptagonal cells pair up to create strand-like connected components (see inset in the rightmost landscape). Defect components of size k are presented in different colours for each k in Figure 2(b), which are captured in the converged configuration at step t = 10 4 shown in the rightmost landscape in (a). Then we compare the number of defect components at the selected iteration steps t = 10, t 0 and 10 4 of the Lloyd process as above (Figure 2(c)). The values in the plot are the mean numbers of defect components of size k ≤ 10 at each of the selected steps, taking 20 realisations into account. At the beginning of the Lloyd process, most defective cells are scattered across the systems as if they repel each other to avoid formation of defective clusters. However, it becomes evident around step t 0 that the number distribution of defect aggregates change into a non-monotone graph in favor of having defects of k = 2 the most, followed by single isolated defects (k = 1) and quadrupole defects (k = 4). The peak at k = 2 is further enhanced as systems are frozen into a defective phase around step t = 10 4 .

Disk packing and hyperuniform systems
We have shown in section 3.1 that the Poisson point pattern undergoes a dynamical arrest beyond the point of topological charge neutrality (t = t 0 ). This dynamical arrest of points is equivalent to the jamming feature in granular systems [18,19]. The hyperuniform 2D landscape represented by n-sided Voronoi cells for n = 5, 6, 7, denoted by C n , can lend itself immediately to the investigation of the 2D disk packing that can be uniquely constructed from this Voronoi landscape. We achieve this by finding the maximum disk of radius r M inscribed in each polygon [13]. A parameter that characterises the cage formed by all the neighbouring disks is the local packing fraction, computed as the disk area (πr 2 M ) divided by the area of the corresponding Voronoi polygon. The harmonic mean of all local packing fractions in the systems corresponds to For t = 10 4 , we consider defective cells denoted by C n (n = 5, 7), highly distorted hexagons C * 6 , which make up the first shell of neighbours to the defective cells, and the remaining hexagons C 6 \ C * 6 . The mean value of packing fraction for each cell type is φ(C 5 ) = 0.791, φ(C 7 ) = 0.805, φ(C * 6 ) = 0.800 and φ(C 6 \ C * 6 ) = 0.845.
the global packing fraction φ. For a perfect 2D hexagonal system, we obtain φ ≈ 0.904, which is acceptably close to its theoretical value of φ lattice ≈ 0.907. We compute the distribution of packing fractions φ at two critical stages of the Lloyd algorithm, namely at the onset of jamming (t = t 0 ) and the final hyperuniform state (t = 10 4 ) (see Figure 3). The results show that the increase in Voronoi cell symmetry enforced by the Lloyd algorithm increases the global packing fraction φ(t) as the system approaches the final hyperuniform state: φ(t 0 ) ≈ 0.809 and φ(10 4 ) ≈ 0.832. Moreover, as the system approaches t = 10 4 , we can observe a skewness of the distribution to the right due to a large number of hexagonal cells C 6 \ C * 6 that are lying out of the first shell of defect components. On the other hand, we can verify that the size polydispersity has been reduced by comparing the standard deviation σ R (t) of the radii distribution of disks: σ R (t 0 ) ≈ 2.028 and σ R (10 4 ) ≈ 1.524.

Discussion
The defective landscape at the end of Lloyd algorithm is a strategic vista for homogenising a point arrangement in 2D: To minimise the total quantizer energy with a fixed number of points, it is inevitable to have pentagons, which are energetically advantageous. But since pentagons alone cannot tile the plane, they pair up with heptagons, each of which area is larger than that of each pentagon as a consequence of sharing an edge. Hence it can be understood as a result of a tradeoff between energy minimisation and geometric frustration in planar systems. At the same time, defect components show a strong preference of forming dipole defects, the simplest possible aggregations, to stabilise the system into a deep local minimum energy scheme. Furthermore, the 2D disk packing constructed from the polygonal tiling of the plane based on the Lloyd algorithm shows that the final defective hyperuniform structure yields a polydisperse packing of disks due to the presence of 5-and 7-coordinated defects (C 5 and C 7 ) as well as distorted hexagons (C * 6 ). This optimal packing of disks obtained by the Lloyd optimisation process has a global packing fraction of φ ≈ 0.832, which is lower than φ ≈ 0.907 for monodisperse disks in a perfect 2D hexagonal lattice. Although it is expected for packing density to increase with polydispersity, this is not observed in our systems. We attribute this result to the Lloyd process, which results in a large number of distorted hexagonal cells having lower local packing density, compared to the regular hexagon found in the 2D lattice system.