Solids are stressed to their breaking point when growing crack networks percolate through the material 1, 2 . Failure by fragmentation may be catastrophic 1, 3 (Fig. 1 ), but this process is also exploited in industrial applications 4 . Moreover, fragmentation of rock and ice is pervasive within planetary shells 1, 5, 6 , and creates granular materials that are literally building blocks for planetary surfaces and rings throughout the solar system [6] [7] [8] [9] [10] (Fig. 1 ). Plato postulated that the idealised form of Earth's building blocks is a cube, the only space-filling Platonic solid 11, 12 . We now know that there is a zoo of geometrically permissible polyhedra associated with fragmentation 13 (Fig. 2) . Nevertheless, observed distributions of fragment mass [14] [15] [16] [17] and shape [18] [19] [20] [21] are self-similar, and models indicate that geometry (size and dimensionality) matters more than energy input or material composition 16, 22, 23 in producing these distributions.
Fragmentation tiles the Earth's surface with telltale mosaics. Jointing in rock masses forms three-dimensional (3D) mosaics of polyhedra, often revealed to the observer by 2D planes at outcrops (Fig. 2) . The shape and size of these polyhedra may be highly regular, even approaching Plato's cube, or resemble a set of random intersecting planes 24 . Alternatively, quasi-2D patterns such as columnar joints sometimes form in solidification of volcanic rocks 25 . These patterns have been reproduced in experiments of mud and corn starch cracks, model 2D fragmentation systems, where the following have been observed: fast drying produces strong tension that drives the formation of primary (global) cracks that criss-cross the sample and make "X" junctions [25] [26] [27] (Fig. 3) ; slow drying allows the formation of secondary cracks that terminate at "T" junctions 26 ; and "T" junctions rearrange into "Y" junctions 25, 28 to either maximise energy release as cracks penetrate the bulk [29] [30] [31] , or during reopening-healing cycles from wetting/drying 32 (Fig.  3) . Whether in rock, ice or soil, the fracture mosaics cut into stressed landscapes (Fig. 3) , form pathways for focused fluid flow, dissolution and erosion that further disintegrate these materials 33 .
Experiments and simulations provide anecdotal evidence that the geometry of fracture mosaics is genetically related to the formative stress field 34 . It is difficult to determine, however, if similarities in fracture patterns among different systems are more than skin deep. First, different communities use different metrics to describe fracture mosaics and fragments, inhibiting comparison among systems and scales. Second, we do not know whether different fracture patterns represent distinct universality classes, or are merely descriptive categories applied to a pattern continuum. Third, it is unclear if and how 2D systems map to 3D.
Here we introduce the mathematical framework of convex mosaics 35 to the fragmentation problem. This approach relies on two key principles: that fragment shape can be well approximated by convex polytopes 24 (2D polygons, 3D polyhedra; Fig. 2a ); and that these shapes must fill space without gaps, since fragments form by the disintegration of solids. Without loss of generality (SI Section 1.1), we choose a model that ignores the local texture of fracture interfaces 36, 37 . Fragments can then be regarded as the cells of a convex mosaic 35 , which Table S2 . may be statistically characterised by three parameters. Cell degree (v) is the average number of vertices of the polytopes, and nodal degree (n) is the average number of polytopes meeting at one vertex 38 : we call [n,v] the symbolic plane. We define the third parameter 0 ≤ p ≡ N R /(N R + N I ) ≤ 1 as the regularity of the mosaic. N R is the number of regular nodes in which cell vertices only coincide with other vertices, corresponding in 2D to "X" and "Y" junctions with n = 4 and 3, respectively. N I is the number of irregular nodes where vertices lie along edges (2D) or faces (3D) of other cells, corresponding in 2D to "T" junctions of n = 2 (Fig. 2b) . We define a regular (irregular) mosaic as having p = 1 (p = 0). For 3D mosaics we also introducef as the average number of faces. In contrast to other descriptions of fracture networks 24 , our framework does not delineate stochastic from deterministic mosaics; networks made from random or periodic fractures may have identical parameter values (Fig. 3) . This theory provides a global chart of geometrically admissible 2D and 3D mosaics in the symbolic plane.
In this paper we measure the geometry of a wide variety of natural 2D fracture mosaics and 3D rock fragments, and find that they form clusters within the global chart. Remarkably, the most significant cluster corresponds to the "Platonic attractor": fragments with cuboid averages. Discrete Element Method (DEM) simulations of fracture mechanics show that cuboid averages emerge from primary fracture under the most generic stress field. Geometric simulations show how secondary fragmentation by binary breakup drives any initial mosaic toward cuboid averages.
The geometric theory of 2D convex mosaics is essentially complete 35 and is given by the formula 38 :
which delineates the admissible domain for convex mosaics within the [n,v] symbolic plane (Fig. 3 ) -i.e., the global chart. Boundaries on the global chart are given by: (i) the p = 1 and p = 0 lines; and (ii) the overall constraints that the minimal degree of regular nodes and cells is 3, while the minimal degree of irregular nodes is 2. We constructed geometric simulations of a range of stochastic and deterministic mosaics (see SI Section 2) to illustrate the continuum of patterns contained within the global chart ( Fig. 3 ) We describe two important types of mosaics, which help to organise natural 2D patterns. First are primitive mosaics, patterns formed by binary dissection of domains. If the dissection is global we have regular primitive mosaics (p = 1) composed entirely of straight lines which, by definition, bisect the entire sample. These mosaics occupy the point [n,v] = [4, 4] in the symbolic plane 35 . In nature, the straight lines appear as primary, global fractures. Next, we consider the situation where the cells of a regular primary mosaic are sequentially bisected locally. Irregular (T-type) nodes are created resulting in a progressive decrease p → 0 and concomitant decreasē n → 2 toward an irregular primitive mosaic. The valuev = 4, however, is unchanged by this process (Fig. 3 ) so in the limit we arrive at [n,v] = [2, 4] . In nature these local bisections correspond to secondary fracturing 3, 34 . Fragments produced from primary vs. secondary fracture are indistinguishable. Further, any initial mosaic subject to secondary splitting of cells will, in the limit, produce fragments withv = 4 (SI Section 1). Thus, we expect primitive mosaics associated with the linev = 4 in the global chart to be an attractor in 2D fragmentation, as noted by 39 , and we expect the average angle to be a rectangle 26 (Fig. 2) . We call this the Platonic attractor. As a useful aside, a planar section of a 3D primitive mosaic (e.g., a rock outcrop) is itself a 2D primitive mosaic (Fig. 2) . The second important pattern is Voronoi mosaics which are, in the averaged sense, hexagonal tilings [n,v] = [3, 6] . They occupy the peak of the 2D global chart (Fig. 3 ).
We measured a variety of natural 2D mosaics (SI Section 2) and found, encouragingly, that they all lie within the global chart permitted by Eq. 1. Mosaics close to the Platonic (v = 4) line include patterns known to arise under primary and/or secondary fracture: jointed rock outcrops, mud cracks, and polygonal frozen ground. Mosaics close to Voronoi include mud cracks and, most intriguingly, Earth's tectonic plates. Hexagonal mosaics are known to arise in the limit for systems subject to repeated cycles of fracturing and healing 25 (Fig. 3) . We thus consider Voronoi mosaics to be a second important attractor in 2D. Horizontal sections of columnar joints also belong to this geometric class; however, their evolution is inherently 3D, as we discuss in section 2.
It is known that Earth's tectonic plates meet almost exclusively at "Y" junctions; there is debate, however, about whether this "Tectonic Mosaic" formed entirely from surface fragmentation, or contains a signature of the structure of mantle dynamics underneath 5, 40, 41 . We examine the tectonic plate configuration 41 as a 2D convex mosaic, treating the Earth's crust as a thin shell. We find [n,v] = [3.0, 5.8], numbers that are remarkably close to a Voronoi mosaic. Indeed, the slight deviation from [n,v] = [3, 6] is because the Earth's surface is a spherical manifold, rather than planar (SI Section 2.3). While this analysis doesn't solve the surface/mantle question, it suggests that the Tectonic Mosaic has evolved from episodes of brittle fracture and healing.
The rest of our observed natural 2D mosaics plot between the Platonic and Voronoi attractors (Fig. 3 ). We suspect that these fractured landscapes, which include mud cracks and permafrost, either: initially formed as regular primitive mosaics and are in various stages of evolution toward the Voronoi attractor; or were Voronoi mosaics that are evolving via secondary fracture towards the Platonic attractor.
There is no formula for 3D convex mosaics analogous to the p = 1 line of Eq. 1 that defines the global chart. There exists a conjecture, however, with a strong mathematical basis 38 ; at present this conjecture extends only to regular mosaics. We define the harmonic degree ash =nv/(n +v). The conjecture is that d <h ≤ 2 d−1 , where d is system dimension. For 2D mosaics we obtain the known result 35, 38h = 2, consistent with the p = 1 substitution in Eq. 1. In 3D the conjecture is equivalent to 3 <h ≤ 4, predicting that all regular 3D convex mosaics live within a narrow band in the symbolic [n,v] plane 38 (Fig. 4) . Plotting a variety of well-studied periodic and random 3D mosaics (SI Section 3), we confirm that all of them are indeed confined to the predicted 3D global chart (Fig. 4) . Unlike the 2D case, we cannot directly measurē n in most natural 3D systems. We can, however, measure the polyhedral cells: the average numbersf ,v of faces and vertices, respectively. Values for [f ,v] may be plotted in what we call the Euler plane, where the lines bounding the permissible domain correspond to simple polyhedra (upper) where vertices are adjacent to three edges and three faces, and their duals which have triangular faces (lower; Fig. 4 ). Simple polyhedra arise as cells of mosaics in which the intersections are generic -i.e., at most three planes intersect at one point -and this does not allow for odd values of v.
As in 2D, 3D regular primitive mosaics are created by intersecting global planes. These mosaics occupy the point [n,v] = [8, 8] on the 3D global chart (Fig. 4) . Cells of regular primitive mosaics have cuboid averages [f ,v] = [6, 8] 35, 38 . This is the Platonic attractor, marked by thev = 8 line in the global chart. 3D Voronoi mosaics, similar to their 2D counterparts, are associated with the Voronoi tessellation defined by some random process. If the latter is a Poisson process then we obtain 35 (Fig. 4) .
Prismatic mosaics are created by regarding the 2D pattern as a base, that is extended in the normal direction. The prismatic mosaic constructed from a 2D primitive mosaic has cuboid averages, and is therefore statistically equivalent to a 3D primitive mosaic. The prismatic mosaic created from a 2D Voronoi base is what we call a columnar mosaic, and it has distinct statistical properties: [n,v] = [6, 12] , [f ,v] = [8, 12] . Thus, the three main natural extensions of the two dominant 2D patterns are 3D primitive, 3D Voronoi, and columnar mosaics.
Regular primitive mosaics appear to be the dominant 3D pattern resulting from primary fracture of brittle materials 42 . Moreover, dynamic brittle fracture produces binary breakup in secondary fragmentation 23, 43 , driving the 3D averages [f ,v] towards the Platonic attractor. The most common example in nature is fractured rock (Figs. 2, 4) . The other two 3D mosaics are more exotic, and seem to require more specialised conditions to form in nature. Columnar joints like the celebrated Giants Causeway, formed by the cooling of large basaltic rock masses 25, 31, 44 (Fig. 4) , appear to correspond to columnar mosaics. In these systems, the hexagonal arrangement and downward (normal) penetration of cracks arise as a consequence of maximising energy release [29] [30] [31] . The only potential example of 3D Voronoi mosaics that we know of is septarian nodules, such as the famous Moeraki Boulders 45 (Fig. 4) . These enigmatic concretions have complex growth and compaction histories, and contain internal cracks that intersect the surface 46 . Similar to primitive mosaics, the intersection of 3D Voronoi mosaics with a surface is a 2D Voronoi mosaic.
We hypothesise that primary fracture patterns are genetically linked to distinct stress fields, in order of most generic to most rare. In a 2D homogeneous stress field we may describe the stress tensor with eigenvalues |σ 1 | ≥ |σ 2 | and characterise the stress state by the dimensionless parameter µ = σ 2 /σ 1 , whose admissible domain is µ ∈ [−1, 1]. In 3D this corresponds to eigenvalues |σ 1 | ≥ |σ 2 | ≥ |σ 3 |, stress state µ 1 = σ 2 /σ 1 , µ 2 = σ 3 /σ 1 , i = sgn(σ 1 ), and domain µ 1 , µ 2 ∈ [−1, 1], |µ 1 | ≥ |µ 2 | (and it is double covered due to i = ±1) (Fig. 5 ). There is a unique map from these 3/8 S10 ). The 2D pattern generator is described by the single scalar functionv(µ) (andn can be computed from [1] ), while the 3D pattern generator is characterised by scalar functions
Computing the full pattern generator is beyond the scope of the current paper, even in 2D. Instead, we perform generic DEM simulations 48, 49 of a range of scenarios to interpret the important primary mosaic patterns described above (see Methods). In 2D, we find that pure shear produces regular primitive mosaics (Fig. 3) , implyingv(−1) ≈ 4. This corresponds to the Platonic attractor. In contrast, hydrostatic tension creates regular Voronoi mosaics ( Fig. 3 ; SI Section 4 ) such that v(1) ≈ 6 -the Voronoi attractor. Both are in agreement with our expectations.
In 3D, we first conducted DEM simulations of hard materials at [µ 1 , µ 2 ] locations corresponding to shear, uniform 2D tension and uniform 3D tension ([−0.5, −0.25], [1, −0.2] and [1, 1]), respectively (Methods, SI Section 4). The resulting mosaics displayed the expected fracture patterns for brittle materials: primitive, columnar and Voronoi, respectively (Fig. 5) . To obtain a global, albeit approximate, picture of the 3D pattern generator we ran additional DEM simulations that uniformly sampled the stress space on a 9 × 9 (∆µ=0.25) grid, for both i= + 1 and i= − 1 (SI Section 4). The constructed pattern generator demarcates the boundaries in stress-state space that separate the three primary fracture patterns (Fig.  5) . The vast proportion of this space is occupied by primitive mosaics, which are also the only pattern generated under negative volumetric stress. Such compressive stress conditions are pervasive in natural rocks. Columnar mosaics are a distant second in terms of frequency of occurrence; they occupy a narrow stripe in the stress space. Most rare are Voronoi mo- (13) Table S2 . Patterns (22) (23) are generated by generic DEM simulation (see Methods): (22) general stress state with eigenvalues σ 1 > σ 2 ; and (23) isotropic stress state with
saics which only occur in a single corner of the stress space (Fig. 5) . Boundaries separating the three patterns shifted somewhat for simulations that used softer materials (Fig. S9 ), but the ranking did not. These primary fracture mosaics serve as initial conditions for secondary fracture. While our DEM simulations do not model secondary fracture, we remind the reader that binary breakup drives any initial mosaic toward an irregular primitive mosaic with cuboid averages (SI Section 1) -emphasising the strength of the Platonic attractor.
Based on the pattern generator (Fig. 5) we expect that natural 3D fragments should have cuboid properties on average, [f ,v] = [6, 8] . To test this we collected 556 particles from the foot of a weathering dolomite rock outcrop (Fig. 6) and measured their values of f and v, plus mass and additional shape descriptors (see Methods; SI Section 5). We find striking agreement: the measured averages [f ,v] = [6.63, 8.93] are within 12 % of the theoretical prediction, and distributions for f and v are centred around the theoretical values. Moreover, odd values for v are much less frequent than even values, illustrating that natural fragments are well approximated by simple polyhedra (Fig. 6) . We regard these results as direct confirmation of the hypothesis, while also recognising significant To better understand the full distributions of fragment shapes, we used geometric simulations of regular and irregular primitive mosaics. The cut model simulates regular primitive mosaics as primary fracture patterns by intersecting an initial cube with global planes (Fig.6 ) while the break model simulates irregular primitive mosaics resulting from secondary fragmentation processes. We fit both of these models to the shape descriptor data using three parameters: one for the cutoff in the mass distribution, and two accounting for uncertainty in experimental protocols (see Methods; SI Section 5). The best fit model, which corresponds to a moderately irregular primitive mosaic, produced topological shape distributions that are very close to those of natural fragments (mean values [f ,v] = [6.58, 8.74]). We also analysed a much larger, previously collected data set (3728 particles) containing a diversity of materials and formative conditions 19 . Although values for v and f were not reported, measured values for classical shape descriptors 19, 21 could be used to fit to the cut and break models (SI Section 5). We find very good agreement (R 2 >0.95), providing further evidence that natural 3D fragments are predominantly formed by binary breakup (SI, Fig.12 ). Finally, we use the cut model to demonstrate how 3D primitive fracture mosaics converge asymptotically toward the Platonic attractor as more fragments are produced (Fig. 6 ). Table S2 .
The application and extension of the theory of convex mosaics provides a new lens to organise all fracture mosaicsand the fragments they produce -into a geometric global chart. There are attractors in this global chart, arising from the mechanics of fragmentation. The Platonic attractor prevails in nature because binary breakup is the most generic fragmentation mechanism, producing averages corresponding to quadrangle cells in 2D and cuboid cells in 3D. Remarkably, a geometric model of random intersecting planes can accurately reproduce the full shape distribution of natural rock fragments. Our findings illustrate the remarkable prescience of Plato's cubic Earth model. One cannot, however, directly 'see' Plato's cubes; rather, their shadows are seen in the statistical averages of many fragments. The relative rarity of other mosaic patterns in nature make them exceptions that prove the rule. Voronoi mosaics are a second important attractor in 2D systems such as mud cracks, where healing of fractures reorganises junctions to form hexagonal cells. Such healing is rare in natural 3D systems. Accordingly, columnar mosaics arise only under specific stress fields, that are consistent with iconic basalt columns experiencing contraction under directional cooling. 3D Voronoi mosaics require very special stress conditions, hydrostatic tension, and may describe rare and poorly understood concretions known as septarian nodules. We have shown that Earth's tectonic mosaic has a geometry that likely arose from brittle fracture and healing, consistent with what is known about plate tectonics 5 (Fig. 3) . This opens the possibility of inferring stress history from observed fracture mosaics. Space missions are accumulating an evergrowing catalogue of 2D and 3D fracture mosaics from diverse planetary bodies, that challenge understanding (Fig. 1) . Geometric analysis of surface mosaics may inform debates on planetery dynamics, such as whether Pluto's polygonal surface (Fig. 1c) is a result of brittle fracture or vigorous convection 7 . Another potential application is using 2D outcrop exposures to estimate the 3D statistics of joint networks in rock masses, which may enhance prediction of rock fall hazards and fluid flow 50 .
Initial samples were randomised cubic assemblies of spheres glued together, with periodic boundary conditions in all directions. The glued contact was realised by a flat elastic cylinder connecting the two particles which was subject to deformation from the relative motion of the glued particles. Forces and torques on the particles were calculated based on the deformation of the gluing cylinder. The connecting cylinder broke permanently if the stress acting upon it exceeded the Tresca criterion 49 . Stress field was implemented by slowly deforming the underlying space. In order to avoid that there is only one percolating crack, we have set a strong viscous friction between the particles and the underlying space. This acts as a homogeneous drag to the particles which ensures a homogeneous stress field in the system.
For any given shear rate the fragment size is controlled by the particle space viscosity and the Tresca criterion limit. We set values that produce reasonable-sized fragments relative to our computational domain, allowing us to characterise the mosaics. Another advantage of the periodic system was that we could avoid any wall effect that would distort the stress field. We note here that it is possible to slowly add a shear component to the isotropic tensile shear test and obtain a structure which has average values ofn andv that are between the primitive mosaics and the Voronoi case. Details of mechanical simulations are discussed in SI Section 4.
field data In the simulation we first computed a regular primitive mosaic by dissecting the unit cube with 50 randomly chosen planes, resulting in 6 × 10 5 fragments. We refer to this simulation as the cut model. Subsequently we further evolved the mosaic by breaking individual fragments. We implemented a standard model of binary breakup 14, 19 to evolve the cube by secondary fragmentation: at each step of the sequence, fragments either break with a probability p b into two pieces, or keep their current size until the end of the process with a probability 1 − p b . The cutting plane is placed in a stochastic manner by taking into account that it is easier to break a fragment in the middle perpendicular to its largest linear extent. Inspired by similar computational models 19 , we used p b dependent on axis ratios (see SI Section 5). This computation, which we call the break model, provides an approximation to an irregular primitive mosaic; this secondary fragmentation process influenced the nodal degreen, but not [v,f ].
In order to compare numerical results with the experimental data obtained by manual measurements, we have to take into account several sampling biases. First, there is always a lower cutoff in size for the experimental samples. We implemented this in simulations by selecting only fragments with m > m 0 , m 0 being the cutoff threshold. Second, there is experimental uncertainty when determining shape descriptors -especially marginally stable or unstable equilibria for the larger dataset (see SI Section 5). We implemented this in the computations by letting the location of the centre of mass be a random variable with variation σ 0 chosen to be small with respect to the smallest diameter of the fragment. We kept only those equilibria which were found in 95% of the cases. Third, there is experimental uncertainty in finding very small faces. We implemented this into the computations by assuming that faces smaller than A 0 P will not be found by experimenters, where P denotes the smallest projected area of the fragment. Using the above three parameters we fitted the seven computational histograms to the seven experimental ones by minimising the largest deviation, and we achieved matches with R 2 max ≥ 0.95 from all histograms (see SI Section 5 for details). Results for the small data set are shown in Fig. 6. 6/8