Hexagulation numbers: the magic numbers of equal spheres on triply periodic minimal surfaces

Regular structures of equal spheres on the triply periodic minimal surfaces known as primitive (P), gyroid (G) and diamond (D) surfaces are enumerated as obtained through Monte Carlo simulations of hard spheres undergoing the Alder transition. Remarkably, there exist magic numbers producing the regular structures, which are simply explained by means of hexagulation numbers defined as H=h2+k2-hk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H=h^2+k^2-hk$$\end{document}, in analogy with the Caspar and Klug’s triangulation numbers, T=h2+k2+hk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=h^2+k^2+hk$$\end{document} for icosahedral viruses, where h and k are equal to nonnegative integers. Understanding the significance of symmetry of the surfaces, the total number of spheres per cubic unit cell N is represented by N=8H\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=8H$$\end{document}, 16H, and 32H for P-, G- and D-surfaces, respectively. Accordingly, these arrangements are analyzed in terms of space groups, equivalent positions (Wyckoff positions), and polygonal-tiling representations. The key is that there is only a limited number of efficient physical design possible even on the triply periodic minimal surfaces.


Introduction
On a flat surface, the hexagonal arrangement is a ubiquitous regular arrangement of spherical particles, arising from dense packing, space division, entropic ordering, or interactions between particles [1]. What is the regular arrangement of particles when the surface is curved? On a spherical surface, this question was firstly raised by Thomson [2], followed by Fejes Tóth [3] and Mackay et al. [4]. For regular arrangements in particular, Goldberg [5] elucidated regular polyhedra consisting of pentagonal and hexagonal faces whose centers correspond to the positions of particles, and later for biological icosahedral viruses, Caspar and Klug proposed a construction principle of regular arrangements inspired by Buckminster Fuller's geodesic dome [6][7][8].
It is well known that spherical virus capsids employ efficient designs using the icosahedral symmetry, leading Dedicated to Professor Alan Mackay on the occasion of his 90th birthday.
to the regularity of the numbers of subunits on the spherical surface, which was explained by Caspar and Klug using the concept of the triangulation number T [6,7]. Surprisingly, we have found only discrete numbers undergoing clear fluid-solid transitions in our MC simulations for the TPMSs. For cubic symmetry in particular, we unravel these magic numbers in terms of a similar number, the hexagulation number H, which is the main result of the present paper. Future applications of such regular arrangements are to construct complex membranes using the concepts of bijels, colloidosomes, and polymersomes [30][31][32][33].
The paper is organized as follows: Section ''Monte Carlo simulations'' briefly explains the Monte Carlo method and presents results of simulations. We list the magic numbers and their space group symmetries and polygonal tilings. All coordinates of regular arrangements are provided in terms of Wyckoff positions of the space groups. In ''Hexagulation Number,'' we first review the triangulation number and then define the hexagulation number to explain the magic numbers on TPMSs. Here, we will focus regular arrangements with original space groups of TPMSs or their cubic subgroups exhibiting a strong similarity with the icosahedral viruses. Finally, discussion is given in the last section.

Monte Carlo simulations
It is widely known that even a purely hard-sphere system undergoes a fluid-solid transition, called the Alder transition, producing the hexagonal arrangement [34], whose entropy-driven ordering mechanism turns out to be universal in soft materials [35]. Above a certain critical density, the solid phase in the hexagonal arrangement gains more entropy than the fluid phase. Do regular arrangements emerge in the solid phase on the TPMSs? Motivated by this possibility, we began an investigation on the G-surface and indeed found the Alder transition, but only for specific numbers of spheres [14]. To provide further insight into the regular arrangement and establish a unified framework, we have extended our MC simulations to the P-and D-surfaces in this study.
The well-known approximate forms [36,37] of the surfaces are described as P : sðqÞ ¼ cos x þ cos y þ cos z ¼ 0; G : sðqÞ ¼ sin x cos y þ sin y cos z þ sin z cos x ¼ 0; D : sðqÞ ¼ sin x sin y sin z þ sin x cos y cos z þ cos x sin y cos z þ cos x cos y sin z ¼ 0; where q ¼ ðx; y; zÞ and the lattice constant is taken to be unity; 2p is omitted in the expressions.
In practice, the centers of spheres are confined between boundaries sðqÞ ¼ À0:1 (green) and sðqÞ ¼ 0:1 (red) as shown in Fig. 1a-c. We assign N a to the number of spheres per cubic unit cell of the surfaces. a stands for P, G, and D. We begin with a random state having a sufficiently small radius of spheres on the TPMSs, perform a MC run without any symmetry input, then increase the radius r with a small increment, and repeat the process. Details of the method of the simulations were given in Ref. [14].
The order parameter h(r) as a function of sphere radius r is defined as where the sum is taken over positions of all sphere centers q i , and hÁ Á Ái implies the MC average. The function f hkl ðqÞ is an invariant function under the operations of the space groups, ex. Im 3m for the P-surface [38]: cosðhxÞ cosðkyÞ cosðlzÞ; where (h, k, l) is a set of integers and C 3 is a group of all cyclic permutations of (x, y, z). Suitable choices of (h, k, l) for each N a give large absolute values of h(r) in the ordered phase. Figure 1d shows a snapshot of an ordered structure for a N P ¼ 72 system on the P-surface, in which alignments are observed in x, y, and z directions despite undulations. To find clear phase behavior, we have employed simulation boxes consisting of multiple cells, here for instance, 3 Â 3 Â 2 cells with periodic boundary conditions; accordingly, the number of spheres is 1296 for N P ¼ 72. The existence of the fluid-solid transition is judged firstly by whether there exists a discontinuous jump or not in the acceptance ratio (AR) curve representing the movable probability of MC trial moves of spheres. If not, the simulation ends up with a frozen random state with a continuous single curve. The transition points are exactly the same as those in Fig. 1e, implying that the ordered phase near the transition region has more movable spaces for spheres than in the disordered phase, and the ordered phase is, therefore, a high entropy phase. We have examined several system sizes to check the results.
Our observation indicates that the number of spheres in a unit cell having a definite fluid-solid transition shows discrete integers. For the P-surface, we find magic numbers N P ¼ 20, 32, 56, 72, and 96; for the G-surface, we find N G ¼ 40, 48, 64, 112, and 144; and for the D-surface, N D ¼ 96, 128, and 224 are obtained. They can be best classified according to the hexagulation numbers H, which will be discussed in the next section. Table 1 summarizes regular structures obtained from MC simulations, where we focus our attention on structures with cubic symmetry, since the TPMSs have inherent cubic space group symmetries, Im 3m, Ia 3d, and Pn 3m for P-, G-and D-surfaces, respectively. Several symmetry operations of the original space groups are broken in some cases, and the ordered arrangements are accordingly assigned to subgroups. For the N P ¼ 72 system, for instance, the arrangement is assigned to the space group Pm 3n & Im 3m. We should mention that the H ¼ 4 class denotes hypothetical structures, while the H ¼ 4 class is obtained from simulations. They are very similar and the difference of the positions on the surfaces is very small, but the H ¼ 4 class is more symmetric than that of H ¼ 4. Because the distances between spheres are more equivalent in H ¼ 4, and consequently in simulations, the symmetry of H ¼ 4 is broken into the entropically more favorable H ¼ 4.
For eager readers, in Table 2 we explicitly provide the values of coordinates of spheres obtained by our MC simulations. One can consult Wyckoff positions in International Tables for Crystallography and evaluate all coordinates. Note that the Wyckoff position 96h of F43c for N D ¼ 96 in the International Tables is not convenient because of the different choices of the origin for symmetry operation. A part of 192h of Fd3c is more useful to reproduce our data, which are presented in Table 3.
Polygonal tilings are constructed from the center positions of spheres, represented in Fig. 2a-i. Since the magic numbers are larger for the G-and D-surfaces, we render five results for the P-surface and two for the G-and Dsurfaces. The set of integers (n 1 ; n 2 ; n 3 ; . . .; n 4 ; n 5 ; n 6 ; . . .;) denotes a tiling of multiple vertex types in the way that n 1gon, n 2 -gon, and n 3 -gon, Á Á Á meet consecutively on each vertex, and superscripts are employed to abbreviate when    Table 3 Struct Chem  International Tables for Crystallography due  possible. A set of integers like (3 6 ; 3 7 ) denotes a tiling composed of two vertex types 3 6 and 3 7 , for instance. The tilings are composed of triangles and quadrangles, where length of all sides is almost equal within a few percent and where quadrangles are not always flat. Unlike the hexagonal arrangement (3 6 ) on a flat plane, it should be noticed that vertices are not equivalent, as is the cases with icosahedral viruses. The first number of the Wyckoff positions listed in Table 1 represents the multiplicity of corresponding vertex types. All these tilings whose Euler characteristic v per cubic unit cell is -4, -8, and -16 for P-, G-, and D-surfaces, respectively, are hyperbolic extensions of the flat plane tessellation (3 6 ) with v ¼ 0. We finally point out strong similarity for the same H classes as shown in Fig. 2d-f (H ¼ 4) and Fig. 2g-i (H ¼ 7).

Hexagulation number
Caspar and Klug showed that icosahedral shells can be constructed from twenty triangles with the triangulation number T adopting the particular integer values 1; 3; 4; 7; 9; 12; . . .; described by with h, k equal to nonnegative integers. See Fig. 3a displaying T-diagram. A lattice point denoted as (h, k) (L ¼ he 1 þ ke 2 ) in the oblique coordinate system spanned by basis vectors e 1 and e 2 on the triangular lattice indicates circled numbers T defined as half the number of lattice points (or equivalently the number of triangles) inside the triangle whose side connects the origin and the point (h, k). Considering the large triangle indicated by the solid line and points occupying the vertices of the triangle in particular, the number of spheres N for the total shell is given by N ¼ 10T þ 2.
By the same token, we define the hexagulation number H as the number of dashed hexagons inside a large hexagon (solid line) on H-diagram as shown in Fig. 3b. Let L ¼ he 1 þ ke 2 be a vector along an edge of the large hexagon, where e 1 and e 2 are basis vectors of the oblique coordinate system. Then, we can easily derive H ¼ jLj 2 . Therefore, we have generating 1; 3; 4; 7; 9; 12; . . .; from sets of integers (h, k). Note that the difference between T and H is just the plus or minus sign of the third term.
According to the space group symmetries, the TPMSs can be constructed from 8, 16, and 32 hexagonal patches per conventional cubic unit cell for P-, G-, and D-surfaces, respectively. Therefore, N a is given by for the P-, G-, or D-surfaces, respectively. The regular structures listed in Table 1 obey these equations. For the special case described by H Ã ¼ 3 in which spheres occupy the vertices of hexagonal patches, we should replace H by H Ã À 1=2 in Eq. (3) as shown in Table 1.
It has been discussed in the literature that hyperbolic geometries on the Poincaré disk can be conformally mapped upon the TPMSs [18,21,39]. See Fig. 4a- On the TPMSs, four hexagonal patches are connected edge to edge around a vertex position of the hexagonal patches making negative curvatures, whose property is rendered on the Poincaré disk (Fig. 4b, d). The radii of circles on the Poincaré disk are constant in each figure, even though they look smaller as it approaches the circumferences. Figure 4a, b displays a system with N P ¼ 56, ðh; kÞ ¼ ð3; 1Þ and H ¼ 7; blue (3 6 ) and red (3 5 :4) circles correspond to 8c and 48i Wyckoff positions of the space group Pn 3n, constructing a polygonal tiling in Fig. 2g. Figure 4c, d illustrates a system with N P ¼ 96, ðh; kÞ ¼ ð4; 2Þ and H ¼ 12. Blue (3 6 ), yellow (3 6 ), and red (3 5 :4) circles correspond to 24g, 24g, and 48h of I 43m, forming a polygonal tiling in Fig. 2c, respectively.

Discussion
We have found that hard spheres on the TPMSs are entropically self-organized when the number of spheres takes some magic numbers N. In analogy with the triangulation number for icosahedral viruses, we have proposed the hexagulation number explaining the magic numbers associated with cubic symmetry. With few exceptions (H Ã ¼ 3 class), N is easily evaluated by H in a unified scheme Eq. (3).
From this scheme, we expect P-G-D triplets, which actually holds for H ¼ 4 and H ¼ 7. However, the first and the second rows of Table 1 appear to be incomplete. In the first row, N P ¼ 20 and N G ¼ 40 are shown, but N D ¼ 80 is missing; any sign of a transition for N D ¼ 80 has not been observed. For N P ¼ 20, the distance between two hard spheres is equidistant (0.3536) and a transition has been found. For N G ¼ 40, the distance is 0.2795 or 0.3062 forming isosceles triangles, whose 8.7 % difference does not matter to undergo a clear transition [14]. However, for N D ¼ 80, the corresponding distance is 0.2165 or 0.2500, whose difference is too large to induce an Alder transition. This incompleteness of the P-G-D triplet indicates that the property of (almost) equidistance between two hard spheres seems to be indispensable to physical ordering even on the curved background.
In the second row of Table 1, we list N G ¼ 48 and N D ¼ 96, but not N P ¼ 24. The physical arrangement is more complex than we expected; surprisingly, for N P ¼ 24 we have had a clear transition, but the unit cell is a doublecell superstructure (N ¼ 192) with rhombohedral symmetry. We have seen a transition to a cubic structure in a simulation box with an odd number of the unit cells in a direction; however, the jump in the acceptance ratio is smaller than that of the transition to the double-cell superstructure, implying that the cubic structure is an artifact by the periodic boundary conditions. This exemplifies the difficulty in simulations, and several box sizes should be tested.
On a sphere, complicated arrangements with lower symmetries have been elucidated [40]. Likewise, our ongoing investigation indeed suggests that there exist complicated regular arrangements with lower symmetries even on the TPMSs. Moreover, not only N P ¼ 24, but also several unprecedented superstructures for other N have been found. They are beyond the scope of the present paper and will be discussed elsewhere.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.