Topological analysis of pattern formation in cooling granular gases confined by elastic wall

For a non-negative integer p, the p-th Betti number of simplicial complexes (see its definition in appendix A) is one of classical topological invariants [20]. It is the number assigned to each simplicial complex, which implies information about its topology. The Betti number implies a 'persistence' with respect to 'continuous deformations'. In other words, the Betti number is invariant under 'continuous deformations'. We note that the Betti number does not imply geometric information (volume, metric etc.) of simplicial complexes. Informally speaking, the p-th Betti number of a simplicial complex X is a number of 'p-dimensional holes' of X. We give the heuristic explanation of the Betti number using figure 1. The simplicial complex in figure 1 consists of '7' vertices and '8' edges.

Zoom In Zoom Out Reset image size

A vertex is called as 0-simplex, and an edge is called as 1-simplex from the simplicial point of view. The simplicial complex in figure 1 has no higher-dimensional simplex such as the face. Now, we compute the 0th Betti number. 0-dimensional hole stands for the connected component following conventions. In short, the 0th Betti number counts the number of connected components. Since the simplicial complex in in figure 1 has one connected component, its 0th Betti number is '1'. On the other hand, a 1-dimensional hole means a circular hole. The simplicial complex in figure 1 has two circular holes along with the hexagon part and the square part. Hence, the 1st Betti number is '2'. These observations can be generalized to the case p ≥ 2. For example, we interpret 2-dimensional hole as a balloon-like hole or 2-sphere-like hole. We introduce the notion of homology in order to deal with such 'p-dimensional holes' in a formal way.

Betti number on the basis of locations of granular particles is calculated using the monotone PNG format with the CHomP [21]. Figure 2 shows the zeroth and first order Betti numbers, which are calculated using monotone PNG. We can readily understand that both B₀ (connection) and B₁ (hole) are equal to '4', as shown in figure 2.

Zoom In Zoom Out Reset image size

For a nonnegative integer p, the p-th Betti number counts 'p-dimensional holes' of a simplicial complex, as we explain in section 2.1. A persistence diagram (PD) is a diagram, which records birth and death of 'p-dimensional holes', i.e. we consider 'time-evolution' of a simplicial complex and are interested in the change of its topology.

The time-evolution, formally speaking, corresponds to filtration of a simplicial complex. A filtration of a simplicial complex X is given by a sequence ${\{{X}_{t}\}}_{t\in {\mathbb{R}}}$ of simplicial sub-complexes ${X}_{t}\subset X$ parametrized by $t\in {\mathbb{R}}$ such that

• ${X}_{{t}_{0}}\subset {X}_{{t}_{1}}$ for t₀ < t₁.
• ${\bigcup }_{t}{X}_{t}=X$ and ${\bigcap }_{t}{X}_{t}=\varnothing$.

A filtration can be understood as a movie film which records how a simplicial complex grows as time evolves. We give an example in figure 3; the filtration ${\{{X}_{t}\}}_{t\in {\mathbb{R}}}$ is an empty complex for sufficiently small t, say t = t₋₁ ; the filtration ${\{{X}_{t}\}}_{t\in {\mathbb{R}}}$ is constant except for t = t₀, t₁, t₂, where it grows up, as shown in figure 3. The p-th PD of a filtration ${\{{X}_{t}\}}_{t\in {\mathbb{R}}}$ is induced by the data of when the p-dimensional holes are born and dead as time evolves. Informally speaking, for a p-dimensional hole h of some X_t, we mark a point $\left({b}_{h},{d}_{h}\right)$ on xy-plane where ${b}_{h}\in {\mathbb{R}}$ (${d}_{h}\in {\mathbb{R}}$, resp.) is the time when the hole h is born (dead, resp.).

Zoom In Zoom Out Reset image size

The p-th PD is obtained by plotting such $\left({b}_{h},{d}_{h}\right)$. The PD is given by multiple set in general since there is a possibility that some holes have the same birth time and death time. We give an example of the PD in figure 4, which is obtained from the filtration in figure 3. Note that the point $\left({t}_{-1},{t}_{0}\right)$ in the PD has multiplicity of '2'. There is an obvious problem with respect to the definition of such a hole h in the time-evolution. The notion of persistent homology and some decomposition theorems are necessary to define PDs in a formal way. Their brief overview is given in appendix C.

Zoom In Zoom Out Reset image size

In this subsection, an overview about how the PD is obtained from point cloud data is given. Since the PD is obtained from a filtration of a simplicial complex, it suffices to construct a filtration of a simplicial complex starting from point cloud data. In practice, a point cloud data is usually given by a finite set S in a Euclidean space ${{\mathbb{R}}}^{N}$. For a positive number ${\mathscr{R}}$, an open covering of S in ${{\mathbb{R}}}^{N}$ given by ${{\mathfrak{U}}}_{{\mathscr{R}}}=\{{{\mathscr{B}}}_{{\mathscr{R}}}(s)\}{}_{s\in S}$ is considered.

Here, ${{\mathscr{B}}}_{{\mathscr{R}}}(s)$ is the ${\mathscr{R}}$-neighborhood of a point $s\in S$ with a radius of ${\mathscr{R}}$. There are several ways to construct a cell complex from the open covering ${{\mathfrak{U}}}_{{\mathscr{R}}}$ : the Čech complex, the Alpha complex, the Vietoris-Rips complex, etc [20, 22, 23].

If we denote one of such complexes by ${C}_{{\mathscr{R}}}$, then, ${\{{C}_{{\mathscr{R}}}\}}_{{\mathscr{R}}\gt 0}$ gives a filtration of a cell complex where we consider ${\mathscr{R}}$ as time, i.e. ${C}_{{{\mathscr{R}}}_{0}}\subset {C}_{{{\mathscr{R}}}_{1}}$ for $0\lt {{\mathscr{R}}}_{0}\lt {{\mathscr{R}}}_{1}$. the PD is obtained from the filtration ${\{{C}_{{\mathscr{R}}}\}}_{{\mathscr{R}}\gt 0}$. We note that the PD is determined by the point cloud data S in ${{\mathbb{R}}}^{N}$.

Figure 5 shows growths of ${{\mathfrak{U}}}_{{\mathscr{R}}}$ of eight point cloud by enlarging ${\mathscr{R}}$. The upper-left frame of figure 5 shows the birth of the hole, when ${\mathscr{R}}={{\mathscr{R}}}_{1}$. The upper-middle frame of figure 5 shows births of eight holes, when ${\mathscr{R}}={{\mathscr{R}}}_{2}$. The upper-right frame of figure 5 shows deaths of holes, when ${\mathscr{R}}={{\mathscr{R}}}_{3}$. Consequently, the PD of these birth-death-sets of holes $({{\mathscr{R}}}_{b},{{\mathscr{R}}}_{d})$ are shown in lower-frame of figure 5, in which ${{\mathscr{R}}}_{b}$ and ${{\mathscr{R}}}_{d}$ correspond to the radius of ${{\mathscr{B}}}_{{\mathscr{R}}}$ (i.e., ${\mathscr{R}}$), which yields the birth and death, respectively. Of course, eight-fold-points are plotted on $({{\mathscr{R}}}_{b},{{\mathscr{R}}}_{d})=({{\mathscr{R}}}_{2},{{\mathscr{R}}}_{3})$ owing to the birth-death of eight holes.

Zoom In Zoom Out Reset image size

The Voronoi's analysis (VA), bond-angle-analysis (BAA) and polyhedral template matching (PTM) have been used to analyze the structure of the crystal on the basis of coordination of atoms. In particular, the analytical tools are freely provided by the codes of Voro++ [24] and Vorotop [25] in Ovito [26]. Then, we can utilize such analyzers in order to analyze the crystal structures of granular particles. Of course, the crystal structure usually postulates the densely packed status of granular particles, where thermal fluctuations of atoms are negligible. We, however, focus on the crystal characteristics, before and after granular particles are densely condensed around the wall owing to their markedly low kinetic energy. Here, schematics of the VA, BAA and PTM are demonstrated, briefly.

In order to discuss the VA, BAA and PTM, we start our discussion by mentioning to the Voronoi's diagram (VD).

${\mathbb{E}}$ is space such as ${\mathbb{E}}\subseteq {{\mathbb{R}}}^{n}$ ($n\in {\mathbb{N}}$). Let's think elements ${{\boldsymbol{g}}}_{1},{{\boldsymbol{g}}}_{2},\ldots {{\boldsymbol{g}}}_{n}\in {\mathbb{E}}$. Here, ${\bigcup }_{k=1}^{n}{{\boldsymbol{g}}}_{k}:= {\mathbb{G}}$.

We define the distance between the point ${\boldsymbol{P}}\in {\mathbb{E}}$ and ${{\boldsymbol{g}}}_{i}$ as

$$\begin{eqnarray*}&&d\left({\boldsymbol{P}},{{\boldsymbol{g}}}_{i}\right),\end{eqnarray*}$$

Then, the Voronoi's domain $R\left({\mathbb{G}};{{\boldsymbol{g}}}_{i}\right)$ is defined by:

$$\begin{eqnarray}&&R\left({\mathbb{G}};{{\boldsymbol{g}}}_{i}\right):= \left\{{\boldsymbol{P}}\in {\mathbb{E}}| d\left({\boldsymbol{P}},{{\boldsymbol{g}}}_{i}\right)\lt d\left({\boldsymbol{P}},{{\boldsymbol{g}}}_{j}\right),\,j\ne i\right\},\end{eqnarray} \tag{ 1 }$$

From equation (1), ${\mathbb{E}}$ is divided into $R({\mathbb{G}},{{\boldsymbol{g}}}_{i})$ ($i\in \left[1,n\right]\cap {\mathbb{N}}$). We call these divided structures of ${\mathbb{E}}$ by $R({\mathbb{G}},{{\boldsymbol{g}}}_{i})$ as the VD.

The VD is described by the Delaunay tetrahedralization in the case of ${\mathbb{E}}\subseteq {{\mathbb{R}}}^{3}$. The significant problem in the VA is classifications of the VD using the typical crystal-coordination such as the face-centered-cubic (FCC), body-centered-cubic (BCC), hexagonal close-packed (HCP), icosahedral (ICO) and hybrid of FCC and HCP (FCC-HCP). The VA by Voro++ applies Weinberg's algorithm [27] to find the most appropriate crystal-coordination among the FCC, BCC, ICO, FCC-HCP and Other from the VD, in which Other indicates the status of coordination of granular particles, from which any crystal structure is not specified with the FCC, BCC, ICO and FCC-HCP.

Of course, there are other methods in order to categorize the crystal structure such as the Common Neighbor Analysis (CNA) [28], PTM [15], BAA [14] other than the VA. The PTM was proposed by Larsen-Schmidt-Schiøtz [15]. In the PTM, the similarity between two diagrams is evaluated by the Root-Mean-Square Deviation (RMSD), which is defined by [15]

$$\begin{eqnarray}&&{{\ell }}_{{\rm{RMSD}}}\left({\boldsymbol{v}},{\boldsymbol{w}}\right)=\mathop{\min }\limits_{s,{\boldsymbol{Q}}}\sqrt{\displaystyle \frac{1}{\aleph }\sum _{i=1}^{\aleph }{\left|\left|s\left[{\vec{v}}_{i}-\bar{{\boldsymbol{v}}}\right]-{\left({\boldsymbol{Q}}{\left[{\vec{w}}_{i}-\bar{{\boldsymbol{w}}}\right]}^{T}\right)}^{T}\right|\right|}_{2}^{2}},\end{eqnarray} \tag{ 2 }$$

where ${\boldsymbol{Q}}$ is the right handed orthogonal matrix , $\bar{{\boldsymbol{v}}}=(1/\aleph ){\sum }_{i=1}^{\aleph }{\vec{v}}_{i},\bar{{\boldsymbol{w}}}=(1/\aleph ){\sum }_{i=1}^{\aleph }{\vec{w}}_{i}$ are barycenter of ${\boldsymbol{v}}$ and ${\boldsymbol{w}}$. s is the optimal scaling of ${\boldsymbol{v}}$. The methods of findings of s and ${\boldsymbol{Q}}$ are proposed by Horn [29] and Theobald [30]. ℵis the number of the neighboring atoms, and $\aleph =6$ for the simple cubic (SC), 12 for the FCC, 12 for the HCP, 12 for the icosahedral (ICO) and 14 for the BCC are used. The detail of algorithm of the PTM is demonstrated by Larsen-Schmidt-Schiøtz [15]. Finally, we find the smallest RMSD among them which are calculated using the SC, FCC, HCP, ICO, and BCC. Additionally, Other corresponds to the state that any structure is not identified even with the SC, FCC, HCP, ICO, and BCC, when the minimum RMSD is larger than its critical value. Actually, ${\boldsymbol{v}}$ is the vector which indicates the vertex-set of the convex hull formed by N-neighboring atoms and ${\boldsymbol{w}}$ is the vector, which indicates the vertex-set of the convex hull of the reference templates, namely, the SC, FCC, HCP, ICO, and BCC (see figure 6 for convex-fulls of the SC, BCC, FCC, HCP and ICO). Hence, we can calculate the distribution of the RMSD (${f}_{{XX}}\left(t,{{\ell }}_{{\rm{RMSD}}}\right),t\in {{\mathbb{R}}}_{+}$: time), when we specify the crystal structure of the atoms using two categories Other and XX (XX := SC, FCC, HCP, ICO and BCC). Finally, we mention to the BAA [14], briefly. The BAA was proposed by Ackland and Jones [14]. The BAA identifies the crystal structure of atom-i with the FCC, HCP, ICO, and BCC from angles (θ_jik) between two bonds (${{\boldsymbol{r}}}_{{ij}}$ and ${{\boldsymbol{r}}}_{{ik}}$), which connect two sets of two neighboring atoms (i-j and i-k). The neighboring atoms are searched by the mean square distance of six neighboring atoms around the atom-i. Once the calculation of the angle between two bonds for several sets of two bonds is finished, functional of χ, which are calculated by angles θ_jik, identifies the crystal structure of atom-i with the FCC, HCP, ICO, BCC or Other.

Zoom In Zoom Out Reset image size

The event-driven (ED) method has been frequently used in order to simulate a large number (N) of granular particles, whose simulation is difficult using the discrete element method (DEM), when the calculation of the $N\left(N-1\right)/2$ paired force-chain is beyond the computational resource. In particular, the ED method is suitable to simulate the dilute granular gases, in which the effects of deformations of contacting granular particles on the collective motion of all the granular particles are negligible.

The algorithm of the ED method is so simple. Provided that the k-th collision occurs at t = t_k (${t}_{k}\in {{\mathbb{R}}}_{+}$), we search for the occurrence-time of k + 1-th collision, namely, t_k+1, which is calculated by

$$\begin{eqnarray}\begin{array}{rcl} & & {t}_{k+1}:= {t}_{k}+\min \left({\rm{\Delta }}{t}_{k}^{{ij}}\right)\,\,\,\left(i,j\in [1,2,\ldots ,N]\right),\\ & & \left|\left|{{\boldsymbol{x}}}_{i}\left({t}_{k}+{\rm{\Delta }}{t}_{k}^{{ij}}\right)-{{\boldsymbol{x}}}_{j}\left({\rm{\Delta }}{t}_{k}^{{ij}}\right)\right|\right|={r}_{i}+{r}_{j},\\ & & {{\boldsymbol{x}}}_{i}\left({t}_{k}+{\rm{\Delta }}{t}_{k}^{{ij}}\right):= {{\boldsymbol{x}}}_{i}\left({t}_{k}\right)+{{\boldsymbol{v}}}_{i}\left({t}_{k}\right){\rm{\Delta }}{t}_{k}^{{ij}},\end{array}\end{eqnarray} \tag{ 3 }$$

where ${{\boldsymbol{x}}}_{i}(t)\in {{\mathbb{R}}}^{d}$ is the location of the granular particle indexed by i at $t,{{\boldsymbol{v}}}_{i}\in {{\mathbb{R}}}^{d}$ is the velocity of the i-th granular particle, r_i is the radius of the i-th granular particle.

Provided the boundary (wall) is considered, the occurrence-time of the collision between the granular particle and wall is calculated in a similar way to equation (3). Afterward, we compare the next collisional time calculated by two granular particles with that calculated by the granular particle and wall and select smaller t_k+1 as next collisional time.

Once t_k+1 is determined, we revise the locations of all the granular particles using ${\boldsymbol{x}}\left({t}_{k+1}\right)={\boldsymbol{x}}\left({t}_{k}\right)\,+{\boldsymbol{v}}({t}_{k}){\rm{\Delta }}{t}_{k}$ (${\rm{\Delta }}{t}_{k}:= {t}_{k+1}-{t}_{k}$). Finally, we revise the velocities of colliding granular particles or velocity of the granular particle colliding with the wall. Provided that the i-th and j-th granular particles collide, the velocities of colliding granular particles change at t = t_k+1 as follows:

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{v}}}_{i}\left({t}_{k+1}\right) & = & {{\boldsymbol{v}}}_{i}\left({t}_{k}\right)-\displaystyle \frac{1+\epsilon }{2}\left({{\boldsymbol{g}}}_{{ij}}\cdot {{\boldsymbol{n}}}_{{ij}}\right){{\boldsymbol{n}}}_{{ij}},\\ {{\boldsymbol{v}}}_{j}\left({t}_{k+1}\right) & = & {{\boldsymbol{v}}}_{j}\left({t}_{k}\right)+\displaystyle \frac{1+\epsilon }{2}\left({{\boldsymbol{g}}}_{{ij}}\cdot {{\boldsymbol{n}}}_{{ij}}\right){{\boldsymbol{n}}}_{{ij}},\\ {{\boldsymbol{g}}}_{{ij}} & = & {{\boldsymbol{v}}}_{i}\left({t}_{k}\right)-{{\boldsymbol{v}}}_{j}\left({t}_{k}\right),\,\,{{\boldsymbol{n}}}_{{ij}}=\displaystyle \frac{{{\boldsymbol{x}}}_{i}\left({t}_{k}\right)-{{\boldsymbol{x}}}_{j}\left({t}_{k}\right)}{\left|\left|{{\boldsymbol{x}}}_{i}\left({t}_{k}\right)-{{\boldsymbol{x}}}_{j}\left({t}_{k}\right)\right|\right|},\end{array}\end{eqnarray} \tag{ 4 }$$

where $\epsilon \in \left[0,1\right]$ is the restitution coefficient, ${{\boldsymbol{g}}}_{{ij}}$ is the relative velocity and ${{\boldsymbol{n}}}_{{ij}}\in {{\rm{\Omega }}}^{2}$ is the relative location-unit-vector.

The calculation of the velocity of the granular particle post-collision with the elastic wall is calculated in a similar way to equation (4) [17].

The restitution coefficient of the granular discs is set as  = 0.85, and the volume fraction of the granular discs is fixed as ϕ = 7.85 × 10⁻². The length of one side of the square box (SQ boundary) is set as L = 300. Three types of the diameter of granular discs (r_d), namely, r_d = 0.15, 0.21 and 0.47 are considered. As a result, total number of granular discs (N) is set as N = 10⁴ in the case of r_d = 0.47, N = 5 × 10⁴ in the case of d = 0.21 and N = 10⁵ in the case of r_d = 0.15 owing to the constant volume fraction (ϕ = 7.85 × 10⁻²). The initial velocities of the granular discs are randomly distributed in the range of ${v}_{x}\in \left[-1,1\right],{v}_{y}\in \left[-1,1\right]$ and the initial positions of the granular discs are randomly distributed inside the elastic wall, namely, the range of $X\in \left[-0.5L,0.5L\right]$ and $Y\in \left[-0.5L,0.5L\right]$. The time-evolution of granular discs is calculated using the event-driven (ED) method [17], because the conventional molecular dynamics requires the vast calculation time, when forces between $N\left(N-2\right)/2$ paired granular discs are calculated.

Figures 7–9 show time-evolutions of the granular discs inside the SQ boundary at N = 10⁴, 5 × 10⁴ and 10⁵, respectively. Figure 7 indicates that some clusters are formed in the vicinity of the wall and the cluster with the maximum size rotates to clockwise direction along the wall. Such a rotation of granular discs along the wall is similar to the emergence of the ordered collective-motion in biological swarm inside the wall [31]. On the other hand, the granular particles tend to move away from the wall, when the heating wall with the constant temperature is used, as reported by Esipov and Pöschel [13]. The numerical results obtained using N = 5 × 10⁴ and 10⁵ do not indicate such a rotation of the cluster with the maximum size, as shown in figures 8 and 9. Indeed, it is commonly observed in numerical results of N = 10⁴, 5 × 10⁴ and 10⁵ that the several small-clusters aggregate toward the larger clusters via their connections during the time-evolution. Of course, the elastic wall is a mathematical toy model to avoid the freeze of the calculation via the excessive accumulation of granular discs with the small kinetic energy in the vicinity of the wall and long range correlation of the granular discs owing to the use of the periodic boundary condition, which is unfavorable to topological analyses, because the distance between two granular discs in the PH must be modified by considering the periodicity at the boundary. The calculations of the PDs are performed using locations of the center of granular discs, namely, $\left({x}_{i},{y}_{i}\right)$ ($i\in \left[1,N\right]\cap {\mathbb{N}}$) and ${\mathscr{R}}$ (:radius of neighborhood (cover) ${{\mathscr{B}}}_{{\mathscr{R}}}\left({s}_{i}\right)$ of point cloud S (${s}_{i}:= \left({x}_{i},{y}_{i}\right)\in S$)). Therefore, we must remind that the connection in d = 0-persistent homology (d: dimension of PH) is not equivalent to the physical connections due to contacting granular discs. Reminding that the domain, which is occupied by neighborhood of i-th point-cloud, whose center is set as $\left({x}_{i},{y}_{i}\right)$, is expressed with ${{\mathscr{B}}}_{{\mathscr{R}}}\left({x}_{i},{y}_{i}\right)$, the connection between ${{\mathscr{B}}}_{{\mathscr{R}}}\left({x}_{i},{y}_{i}\right)$ and ${{\mathscr{B}}}_{{\mathscr{R}}}\left({x}_{j},{y}_{j}\right)$ is defined by $\left|{{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{j}\right|\lt 2{\mathscr{R}}$ (${{\boldsymbol{x}}}_{i}:= \left({x}_{i},{y}_{i}\right),2{\mathscr{R}}\gt d$). ${\mathscr{R}}$ is increased from zero to $\infty$, continuously in ${\mathscr{R}}\in {{\mathbb{R}}}_{+}$. ${\mathscr{R}}={{\mathscr{R}}}_{b}$ is called as the birth of the hole, when the hole emerges when ${\mathscr{R}}={{\mathscr{R}}}_{b}$, whereas ${\mathscr{R}}={{\mathscr{R}}}_{d}$ is called as the death of the hole, when the hole disappears, when ${\mathscr{R}}={{\mathscr{R}}}_{d}$. The plot of $\left({{\mathscr{R}}}_{b},{{\mathscr{R}}}_{d}\right)$ is called as the PD, as discussed in section 2.2. Additionally, ${\ell }:= {{\mathscr{R}}}_{d}-{{\mathscr{R}}}_{b}$ is called as the life-span of the hole, and defined by:

$$\begin{eqnarray}&&{\ell }:= {{\mathscr{R}}}_{d}-{{\mathscr{R}}}_{b},\end{eqnarray} \tag{ 5 }$$

The image of the life-span (ℓ) is shown in the lower-frame of figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The open source program homcloud [32] is used to calculate the PD in our study. The PD is plotted using not $\left({{\mathscr{R}}}_{b},{{\mathscr{R}}}_{d}\right)$ but $\left({{\mathscr{R}}}_{b}^{2},{{\mathscr{R}}}_{d}^{2}\right)$ owing to the specification of the homcloud. Figure 10 shows PDs (left-half frames) and birth of holes with ${{\ell }}_{\max }$ (maximum value of ℓ) via connections of some of ${{\mathscr{B}}}_{{\mathscr{R}}}\left({s}_{i}\right)$ ($i\in [1,N]$) at N_c = 0, 2.15M, 3.85M and 11.4M (right-half frames), when N = 10⁵ (N_c: collision number, M: million), respectively. The color of the contour expresses the number of holes, which has $\left({{\mathscr{R}}}_{b}^{2},{{\mathscr{R}}}_{d}^{2}\right)$, where $U\left({{\mathscr{R}}}_{b}^{2},{{\mathscr{R}}}_{d}^{2}\right)$ is the square domain, whose center of gravity is $\left({{\mathscr{R}}}_{b}^{2},{{\mathscr{R}}}_{d}^{2}\right)$ and length of one side of the square is set as ${\rm{\Delta }}\in {{\mathbb{R}}}_{+}$. Readers remind that the number of holes which yield ${{\ell }}_{\max }$ is unity in all the frames in figure 10. The magnitude of ${{\mathscr{R}}}_{b}$ obtained using ${{\ell }}_{\max }$ increases in accordance with the increase in N_c, as shown in the right-half frame of figure 10, because the growth of clusters of granular discs tends to enlarge the vacant space. The PD has the clear structure at N_c = 0, whereas such a clear structure in the PD becomes blurred at N_c = 14.45M.

Zoom In Zoom Out Reset image size

The expeditious method to read the characteristics of the PD is to calculate the life-span-distribution (lsd) by each time step. Then, we approximate the lsd, namely, $f\left({\ell }\right):= \left(1/N\right){\sum }_{i=1}^{N}\delta \left({\ell }-{{\ell }}_{i}\right)$ with

$$\begin{eqnarray}&&{f}_{{\rm{ap}}}\left({\ell }\right):= {\mathfrak{A}}\exp \left[-{\mathfrak{B}}{\left({\ell }-{\mathfrak{C}}\right)}^{n}\right]{\left({\mathfrak{D}}{\ell }+{\mathfrak{E}}\right)}^{m}\end{eqnarray} \tag{ 6 }$$

where ${\mathfrak{A}},{\mathfrak{B}}$, and ${\mathfrak{D}}\subseteq {{\mathbb{R}}}_{+},{\mathfrak{C}},{\mathfrak{E}},n$ and m $\subseteq {\mathbb{R}}$.

The reason why we use equation (6) is that it can express the hybrid of the logarithmic convex (log-convex: for $\theta \in [0,1],\theta \mathrm{log}f\left(x\right)+\left(1-\theta \right)\mathrm{log}f\left(y\right)\leqslant \mathrm{log}f\left(\theta x+(1-\theta )y\right)$ ($a\in {{\mathbb{R}}}_{+}$ and x < y)) and logarithmic concave (log-concave: for $\theta \in [0,1],\theta \mathrm{log}f\left(x\right)+\left(1-\theta \right)\mathrm{log}f\left(y\right)\geqslant \mathrm{log}f\left(\theta x+(1-\theta )y\right)$ ($a\in {{\mathbb{R}}}_{+}$ and x < y)) using appropriate set of $\left(n,m\right)$ in equation (6). Similarly, equation (6) is also able to express completely concave or convex for all the range of ℓ (i.e., ${\ell }\in {{\mathbb{R}}}_{+}$).

Figures 11–16 show $f\left({\ell }\right)$ versus ℓ together with ${f}_{{\rm{ap}}}\left({\ell }\right)$ for d = 0 (connection) and 1 (hole) in cases of N = 10⁴, 5 × 10⁴ and 10⁵, whose ${\mathfrak{A}},{\mathfrak{B}},{\mathfrak{C}},{\mathfrak{D}},{\mathfrak{E}},{\mathfrak{F}},n$ and m in equation (6) are defined in table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  ${\mathfrak{A}},{\mathfrak{B}},{\mathfrak{C}},{\mathfrak{D}},{\mathfrak{E}},n$ and m in equation (6) and corresponding figure number, when d = 0 and 1, N = 10⁴, 5 × 10⁴ and 10⁵, and SQ boundary are used.

Form	d	N	Nc	${\mathfrak{A}}$	${\mathfrak{B}}$	${\mathfrak{C}}$	${\mathfrak{D}}$	${\mathfrak{E}}$	n	m	Figure No.
SQ	0	104	0	0.8	0.29	0	1	1	2	0.2	Figure 11
SQ	0	104	2.5M	1.2	0.29	0	1	1	1.8	0.35	Figure 11
SQ	0	104	5M	2.25	0.15	0	1	1	0.8	−2.5	Figure 11
SQ	0	104	12.5M	0.25	0.05	0	1	0.1	1.26	−1.5	Figure 11
SQ	0	104	19.9M	0.15	0.03	0	1	1	1.26	−1.5	Figure 11
SQ	0	5 × 104	0	4.5	5.5	0	1	1	2.2	−2	Figure 13
SQ	0	5 × 104	1M	4.5	5.5	0	1	1	2.2	−2	Figure 13
SQ	0	5 × 104	1.5M	8	3.25	0	1	1	1	−3	Figure 13
SQ	0	5 × 104	8.85M	12	3.25	0	1	1	0.25	−3	Figure 13
SQ	0	105	0	10	25	0	1	1	2	−3	Figure 15
SQ	0	105	2.15M	10	15	0	1	1	1.7	−3	Figure 15
SQ	0	105	3.85M	8	4.2	0	1	0.75	0.5	−5	Figure 15
SQ	0	105	8.85M	0.5	0.1	0	1	0.75	1.2	−5	Figure 15
SQ	0	105	14.5M	0.5	0.1	0	1	0.75	1.2	−5	Figure 15
SQ	1	104	0	0.003	2.5 × 10−5	0	0.067	10−5	4	−1.5	Figure 12
SQ	1	104	1M	0.005	6.1 × 10−9	0	0.2	10−5	2	−1.5	Figure 12
SQ	1	104	5M	0.005	2.2 × 10−10	0	0.2	10−5	3.2	−1.75	Figure 12
SQ	1	104	12.5M	0.005	3.57 × 10−22	0	0.2	10−5	6.2	−1.35	Figure 12
SQ	1	5 × 104	0	5 × 10−3	3.14 × 10−2	0	0.2	10−5	3	−1.5	Figure 14
SQ	1	5 × 104	1M	5 × 10−3	1.6 × 10−2	0	0.2	10−5	3	−1.5	Figure 14
SQ	1	5 × 104	1.5M	5 × 10−3	0.136	0	0.2	10−5	1.5	−1.75	Figure 14
SQ	1	5 × 104	7.4M	1.2 × 10−3	3.79 × 10−3	0	0.2	10−5	0.98	−2.4	Figure 14
SQ	1	105	0	5 × 10−3	0.25	0	0.5	0	3	−1.5	Figure 16
SQ	1	105	2.15M	0.02	1	0	0.4	0	1	−1.5	Figure 16
SQ	1	105	3.85M	0.01	0.71	0	0.4	0	0.5	−2	Figure 16
SQ	1	105	14.5M	0.01	1	0	0.4	0	1.5 × 10−3	−2.5	Figure 16

The authors will consider that one standard for the evaluation of the topological phase-transition of the granular discs is determined by the switch between the log-convex and log-concave of $f\left({\ell }\right)$ (or ${f}_{{\rm{ap}}}\left({\ell }\right)$), as discussed later. Figure 11 indicates that $f\left({\ell }\right)$ (or ${f}_{{\rm{ap}}}\left({\ell }\right)$) follows the log-concave at N_c = 0 and 0.25M, whereas $f\left({\ell }\right)$ (or ${f}_{{\rm{ap}}}\left({\ell }\right)$) follows the log-convex at 0.5M≤ N_c. Consequently, we conclude that the topological phase-transition occurs in the range of 0.25M< N_c < 0.5M, when d = 0 and N = 10⁴. Similarly, figure 12 indicates that $f\left({\ell }\right)$ (or ${f}_{{\rm{ap}}}\left({\ell }\right)$) follows the log-convex at ${\ell }\leqslant 9$ and log-concave at $9\lt {\ell }$ in the case of N_c = 0, whereas $f\left({\ell }\right)$ (or ${f}_{{\rm{ap}}}\left({\ell }\right)$) follows the log-convex at 0 ≤ ℓ. As a result, the topological phase-transition occurs in the range of N_c ∈ [0, 1M], when d = 1 and N = 10⁴. These topological phase-transitions are obtained at N_c ∈ [1M,1.5M] in both cases of d = 0 and 1, when N = 5 × 10⁴, as shown in figures 13 and 14, whereas it is obtained at N_c ∈ [2.15M,3.85M] in both cases of d = 0 and 1, when N = 10⁵, as shown in figures 15 and 16.

The lsd obtained using d = 1-PD, which follows the log-convex, indicates that the number of holes with the large diameter ${{\mathscr{R}}}_{d}$ and small ${{\mathscr{R}}}_{b}$ (i.e., long life-span) is larger than that obtained using the d = 1-PD, which follows the log-concave. The holes with the long life-span are generated by the increase of vacant space and connections of point-clouds by small ${{\mathscr{R}}}_{b}$, which are caused by the growths of clusters. Consequently, the switch between the log-concave and log-convex becomes the standard to measure the drastic growth of clusters, by which the topological characteristics of granular gases changes, markedly.

From above discussions, we confirmed that the topological phase-transition can be judged by the switch of the life-span-distribution (lsd) (i.e., $f\left({\ell }\right)$) between the log-convex and log-concave. Here, we investigate the topological phase-transition using the Betti number, B₀ (the zeroth order) and B₁ (the first order). As discussed in section 2, the Betti number is calculated using the PNG versions of figures 7–9. Therefore, the accuracies of the calculations of B₀ and B₁ depend on the number of pixels in PNG versions of figures 7–9. Figure 17 shows B₀/N and B₁/N in cases of N = 10⁴ (upper-left frame), 5 × 10⁴ (upper-right frame) and N = 10⁵ (lower-left frame). B₀/N ≃ 0.96 at N_c = 0 obtained using N = 10⁵ indicates that the numerical error exists owing to the insufficient number of pixels, whereas B₀/N = 1 at N_c = 0 obtained using N = 10⁴ and 5 × 10⁴ indicate that the number of pixels is adequate to calculate B₀ and B₁ with the good accuracy.

Zoom In Zoom Out Reset image size

Firstly, we confirm that there exist typical three phases (i.e., phases-I, II and III) in cases of N = 5 × 10⁴ and 10⁵. In the phase-I, B₀ (B₁) decreases (increases), drastically, as shown in N_c ∈ [0,0.5M] when N = 5 × 10⁴, and N_c ∈ [0,1M] when N = 10⁵. In the phase-II, B₀ (B₁) changes, slightly, as shown in N_c ∈ [0.5M,1.5M] when N = 5 × 10⁴ and N_c ∈ [1M,2M] when N = 10⁵. Finally, in the phase-III, B₀ (B₁) decreases (increases), drastically, as shown in N_c ∈ [1.5M, $\infty$) when N = 5 × 10⁴ and N_c ∈ [2M, $\infty$) when N = 10⁵. We are, however, unable to identify phases-I and II, clearly, when N = 10⁴, whereas the phase-III is conformed in N_c ∈ [0.25M, $\infty )$, when N = 10⁴. The increase and decrease in B₀ and B₁ in the phase-III are caused by temporal changes in form of rotating cluster along the wall obtained using N = 10⁴, so that such changes in phase-III do not attribute to intrinsic changes in the topological phase, which characterize the pattern formation from the HCS. It is the significant result that the topological phase-transition from the phase-I to the phase-II is identified by not the lsd but B₀ and B₁, whereas the phase-transition due to the topological phase-transition from the phase-II to the phase-III is identified by both the lsd and B₀ and B₁.

Now, we investigate the PH in the case of the CI boundary in order to investigate the effects of form of the elastic wall on the PH. The radius of the granular disc (r_d) is set as r_d = 0.15 and radius of the CI (R_d) is set as R_d = 150. The set of the restitution coefficient and initial velocities of the granular discs are same as those used in calculations in the case of the SQ boundary. The total number of the granular discs (N) is set as $N={10}^{5}$. As a result, the volume fraction of granular discs (ϕ) is calculated as ϕ = 0.1, because all the granular discs are randomly distributed inside the CI boundary at t = 0.

Figure 18 shows snap-shots of the granular discs confined by the CI boundary at N_c = 0 (upper-left frame), N_c = 2.5M (upper-center frame), N_c = 5M (upper-right frame), N_c = 7.5M (lower-left frame), N_c = 15M (lower-center frame) and N_c = 35M (lower-right frame). We can confirm the clear pattern formation at 5M≤N_c. In particular, the clusters seem to grow from the wall at N_c = 35M.

Zoom In Zoom Out Reset image size

Figure 19 shows the d = 1-PDs at N_c = 0 (top-left frame), N_c = 2.5M (top-right frame), N_c = 5M (middle-left frame), N_c = 7.5M (middle-right frame), N_c = 15M (bottom-left frame) and N_c = 35M (bottom-right frame), when N_c = 10⁵ and CI boundary. The tendency of changes in the d = 1-PD during the time-evolution is similar to that obtained using the SQ boundary in figure 10. Figures 20 and 21 show $f\left({\ell }\right)$ versus ℓ together with ${f}_{{\rm{ap}}}\left({\ell }\right)$ (in equation (6)) versus ℓ when d = 0 and d = 1, respectively, where ${\mathfrak{A}},{\mathfrak{B}},{\mathfrak{C}},{\mathfrak{D}},{\mathfrak{E}},{\mathfrak{F}},n$ and m in equation (6) are defined in table 2. We can confirm that the topological phase-transition occurs in N_c ∈ [2.5M,5M] in both cases of d = 0 and 1, as shown in figures 20 and 21.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  ${\mathfrak{A}},{\mathfrak{B}},{\mathfrak{C}},{\mathfrak{D}},{\mathfrak{E}},n$ and m in equation (6) and corresponding figure number, when d = 0 and 1, N = 10⁵ and CI boundary are used.

Form	d	N	Nc	${\mathfrak{A}}$	${\mathfrak{B}}$	${\mathfrak{C}}$	${\mathfrak{D}}$	${\mathfrak{E}}$	n	m	figure No.
CI	0	105	0	1	34.8	0	1	3.15	1.74	2	Figure 20
CI	0	105	2.5M	900	18	0	1	−10−4	0.8	1	Figure 20
CI	0	105	5M	100	9.5	0	1	−10−4	0.48	0.35	Figure 20
CI	0	105	7.5M	50	8	0	1	−10−4	0.4	0.2	Figure 20
CI	0	105	35M	50	8	0	1	−10−4	0.3	0.15	Figure 20
CI	1	105	0	0.05	0.75	0	1	10−7	3	−1.4	Figure 21
CI	1	105	2.5M	0.05	0.75	0	1	10−7	2	−1.4	Figure 21
CI	1	105	5M	0.02	0.15	0	1	10−7	1.2	−1.75	Figure 21
CI	1	105	35M	0.02	0.15	0	1	10−7	1.2	−1.75	Figure 21

Figure 22 shows B₀/N, B₁/N and $\left({B}_{0}+{B}_{1}\right)/N$. Similarly to B₀/N obtained using N = 10⁵ and the SQ boundary, B₀/N ≃ 0.96 indicates that the numerical error exists owing to the insufficient number of pixels. We can readily confirm that phases-II and III exist in N_c ∈ [0,2M] and N_c ∈ [2M, $\infty$). It is not obvious in the present study whether the lack of the phase-I in the case of the CI boundary is caused by the different form of boundary from the SQ boundary or larger volume fraction (ϕ = 0.1) in CI boundary than that (ϕ = 0.078 5) in the SQ boundary. The phase-transition due to the topological phase-transition from the phase-II to the phase-III is obtained using both the lsd and Betti number in the case of the CI boundary.

Next, we investigate the PH in the case of the SP boundary. As discussed in section 2, the 3D physical space postulates the PH with d = 2, which corresponds to the spherical surface (void). Here, we consider two types of N. In one case, N = 5 × 10⁴ granular spheres with r_d = 1.5 are calculated. In the other case, N = 10⁴ granular spheres with r_d = 1.5 are calculated. In both cases, the radius of the SP boundary is set as R_d = 150 and restitution coefficient is set as 0.85. Consequently, the volume fraction is 5% in the case of N = 5 × 10⁴ and 1% in the case of N = 10⁴, because all the granular spheres are randomly (almost homogeneously) distributed inside the SP boundary. The restitution coefficient of the granular spheres is set as 0.85 and the initial velocities of granular spheres are randomly selected in the ranges of ${v}_{x}\in [-1,1],{v}_{y}\in [-1,1]$ and ${v}_{z}\in [-1,1]$.

Zoom In Zoom Out Reset image size

Figure 23 shows the snapshots of the granular spheres confined by the SP boundary, when N = 10⁴ (upper two raw) and N = 5 × 10⁴ (lower two raw). The granular spheres cluster in the vicinity of the SP boundary with the band-like form as N_c increases, when N = 10⁴,. On the other hand, some band-like distributions of granular spheres emerge along the SP boundary in accordance with the increase in N_c when N = 5 × 10⁴. The aggregation of granular spheres in the vicinity of the wall at large N_c is similar to those in 2D cases of the SQ and CI boundaries.

Zoom In Zoom Out Reset image size

Figures 24–27 show time-evolutions of d = 1-PDs and d = 2-PDs obtained using N_c = 10⁴ and 5 × 10⁴, respectively. The tendency of temporal changes in d = 1-PDs is similar to those obtained using the SQ and CI boundaries. d = 2-PDs obtained using N = 10⁴ indicate that the life-spans of most of granular spheres approach to zero at N_c = 5.5M, when N = 10⁴, although the life-spans of some of granular spheres are finite, as observed in d = 1-PD at N_c = 5.5M in figure 24. Meanwhile, the life-spans of granular spheres observed in d = 2-PDs are finite at N_c = 40M, when N = 5 × 10⁴, as shown in figure 27.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 28–33 show $f\left({\ell }\right)$ versus ℓ together with ${f}_{{\rm{ap}}}\left({\ell }\right)$ versus ℓ obtained using N = 10⁴ and 5 × 10⁴, when d = 0, 1 and 2, respectively, where ${\mathfrak{A}},{\mathfrak{B}},{\mathfrak{C}},{\mathfrak{D}},{\mathfrak{E}},{\mathfrak{F}},n$ and m in equation (6) are defined in table 3. Figures 28 and 30 indicate that the topological phase-transition occurs in N_c ∈ [2M,2.5M] in the case of N = 10⁴, whereas all of $f\left({\ell }\right)$ for d = 2 follow log-convex at N_c ∈ [0, $\infty$), when N = 10⁴, as shown in figure 32. Similarly, the topological phase-transition occurs in N_c ∈[3M,4M], when N = 5 × 10⁴, as shown in figure 29. The log-concave of f(ℓ) remains at N_c = 45M, when d = 1 and N = 5 × 10⁴, as shown in figure 31. Then, we cannot determine the topological phase-transition from the lsd for d = 1, when N = 5 × 10⁴. Similarly to f(ℓ) obtained using d = 2 and N = 10⁴, all of $f\left({\ell }\right)$ follow the log-convex in ${N}_{c}\in [0,\infty )$, as shown in figure 33. Consequently, we conclude that the lsd for d = 2 can not be used as the standard, which determines the topological phase-transition in accordance with the switch of the lsd between the log-concave and log-convex in both cases of N = 10⁴ and 5 × 10⁵.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  ${\mathfrak{A}},{\mathfrak{B}},{\mathfrak{C}},{\mathfrak{D}},{\mathfrak{E}},n$ and m in equation (6) and corresponding number of the figure, when d = 0, 1 and 2, N = 10⁴ and 5 × 10⁴ and SP boundary are used.

Form	d	N	Nc	${\mathfrak{A}}$	${\mathfrak{B}}$	${\mathfrak{C}}$	${\mathfrak{D}}$	${\mathfrak{E}}$	n	m	figure No.
SP	0	104	0	0.012	0.003	1	0.8	1	2	0.75	Figure 28
SP	0	104	2M	0.012	0.003	1	0.8	1	2	0.75	Figure 28
SP	0	104	2.5M	0.02	0.003	1	0.8	1	0.95	1	Figure 28
SP	0	104	3M	5 × 107	18	−1	1	0	2	0.15	Figure 28
SP	0	104	5.5M	1	0	0	1	−614	1	0.77	Figure 28
SP	1	104	0	1	2.25 × 10−4	0	3.6	1.83	2	−1.15	Figure 30
SP	1	104	2M	1	$2.25\times {10}^{-4}$	0	3.6	1.83	2	−1.15	Figure 30
SP	1	104	2.5M	1	0.015	0	2.46	1.98	1	−1.32	Figure 30
SP	1	104	3M	1	0.012	0	2.33	2.26	1	−1.45	Figure 30
SP	1	104	3.5M	1	0.1	0	2.13	2.75	1	−1.77	Figure 30
SP	1	104	5.5M	1	0	0	107	7.88	1	−1.08	Figure 30
SP	2	104	0	1	0	0	1.16	2.9	1	−1.95	Figure 32
SP	2	104	2M	1	0	0	1.16	2.9	1	−1.95	Figure 32
SP	2	104	2.5M	1	0	0	1.16	2.9	1	−2	Figure 32
SP	2	104	3M	1	0	0	1.16	2.9	1	−2.2	Figure 32
SP	2	104	5.5M	1	0	0	9.33	4.64	1	−1.1	Figure 32
SP	0	5 × 104	0	1	1.62 × 10−9	−48	1	0	5.5	2	Figure 29
SP	0	5 × 104	2M	1	2.1 × 10−8	−42.3	1	0	5	2	Figure 29
SP	0	5 × 104	3M	1	1.66 × 10−2	−12.3	1	0	2	2	Figure 29
SP	0	5 × 104	4M	1	8.8 × 10−3	-6,15	1	0	2	0.6	Figure 29
SP	0	5 × 104	5M	1	4.2 × 10−3	−0.144	1	0	2.5	−1.8	Figure 29
SP	0	5 × 104	45M	1	2.44 × 10−2	−2.13	1	0	1	−1.5	Figure 29
SP	1	5 × 104	0	0.25	10−7	−3	1	0	5	−1.5	Figure 31
SP	1	5 × 104	2M	0.25	10−7	−3	1	0	5	−1.5	Figure 31
SP	1	5 × 104	3M	0.25	1.6 × 10−6	−3	1	0	4	−1.5	Figure 31
SP	1	5 × 104	4M	0.25	2.5 × 10−6	−3	1	0	3.75	−1.5	Figure 31
SP	1	5 × 104	5M	0.25	2.5 × 10−6	−1	1	0	3.5	−1.75	Figure 31
SP	1	5 × 104	45M	0.25	5 × 10−6	−1	1	0	2.75	−1.5	Figure 31
SP	2	5 × 104	0	0.25	10−5	−3	1	0	3	−2	Figure 33
SP	2	5 × 104	4M	0.25	10−5	−3	1	0	3	−2	Figure 33
SP	2	5 × 104	10M	0.1	10−5	−3	1	0	3	−2	Figure 33
SP	2	5 × 104	45M	0.1	10−5	−3	1	0	3	−2	Figure 33

The Betti number is not calculated in the case of the SP boundary, because the monotone (0,1) pixels, which are used in the calculation of B₀ and B₁ in the SQ and CI boundaries, are unable to consider the depth-direction in 3D.

Figure 34 shows snapshots of the categorized coordination per a granular particle and its fraction obtained using the BAA at N_c = 0, 4M, 5M, 6M, 10M, 20M, 30M, 40M and 47.4M. We can confirm that the spatially random locations of granular spheres at N_c = 0 (t = 0) does not obtain any crystal structure. In short, all the coordination of granular particles is categorized as Other. Meanwhile, the fraction of granular spheres with the coordination categorized as the HCP, FCC and BCC (coordination) increases, as N_c increases. Actually, 20% of granular spheres are categorized as the structured coordination, namely, FCC, HCP, BCC and ICO at N_c = 47.7M. Of course, most of granular spheres condensate around the wall at N_c = 47.7M, so that such a status is not appropriate to be called as the granular gases.

Zoom In Zoom Out Reset image size

Figure 35 shows the fractions of the categories of the coordination versus N_c obtained using the BAA (Other: left frame, FCC, HCP, BCC and ICO: right frame). The enlarged figure (log-log plot of the fraction of Other versus N_c) is added to the left frame of figure 35 in order to show that the fraction of Other decreases by the inverse power law function of N_c. As shown in figure 35, the fraction of the category of the coordination is invariant during N_c < 4M and decreases in accordance with $130/{\left({N}_{c}/M+4\right)}^{-0.125}$. Then, we can discriminate the topological phase of the granular spheres between N_c < 4M (HCS) and 4M≤N_c (growth of clusters). Reminding that the lsd for d = 0 changes from the log-concave to the log-convex at 3M< N_c < 4M, as shown in figure 29, the BAA also specifies the topological phase-transition of the granular spheres with a similar accuracy to the lsd obtained using d = 0-PD.

Zoom In Zoom Out Reset image size

Figure 36 shows the snapshots of the categorized coordination per a granular sphere and its fraction at N_c = 0, 2M 4M, 5M, 6M, 10M, 20M, 30M, 40M and 47.4M obtained using the VA. We can confirm that 3.1% of randomly located granular spheres at N_c = 0 can be categorized as HCP or FCC (coordination). Finally, 32.8% of granular spheres obtain the crystal structures, namely, FCC, HCP, BCC, FCC-HCP and ICO at N_c = 47.7M.

Zoom In Zoom Out Reset image size

Figure 37 shows the fractions of the categories of the coordination versus N_c obtained using the VA (Other: left frame, FCC, BCC, ICO, FCC-HCP and HCP: right frame). The enlarged figure (log-log plot of the fraction of Other versus N_c) is added to the left frame of figure 34 in order to show that the fraction of Other decreases by the inverse power law function of N_c. As shown in figure 34, the fraction of the category Other is almost invariant during N_c < 4M and decreases in accordance with $145/{\left({N}_{c}/M+4\right)}^{-0.2}$. Then, the VA-results indicate that the topological phase-transition of the granular spheres occurs between N_c < 4M (HCS) and 4M≤N_c (growth of clusters) similarly to the PD and BAA. The growths of the fractions of the FCC-HCP and HCP are markedly larger than those of the BCC, ICO and FCC. Then, we can conclude that the HCP type coordination is dominant in both cases of the BAA and VA, when the structured coordination are assigned to granular spheres.

Zoom In Zoom Out Reset image size

We consider the categorization of the coordination obtained using the PTM. Figures 38–46 show granular particles inside the SP boundary, which are colored in accordance with the assigned category of the coordination (i.e., All:=Other, FCC, HCP, BCC, ICO and SC), at N_c = 0, 4M, 5M, 6M, 10M, 20M, 30M, 40M, and 47.7M together with ${\sum }_{{XX}}{f}_{{XX}}\left({{\ell }}_{{\rm{RMSD}}}\right)$ versus ℓ_RMSD (XX:=FCC, HCP, BCC, ICO and SC) in their top-left frames. Additionally, other remained frames except for top-left frame in figures 38–46 show the granular particles inside the SP boundary, which are colored by only two categories Other and XX (XX:=FCC, HCP, BCC, ICO and SC) together with the distribution of the RMSD, namely, ${f}_{{XX}}\left({{\ell }}_{{\rm{RMSD}}}\right)$ versus ${{\ell }}_{{\rm{RMSD}}}$, when all the granular particles are categorized by Other and XX (i.e., FCC, HCP, BCC, ICO and SC). Figure 38 shows that 48.7% of N is categorized as Other, whereas 45.3% of N is categorized as the SC and 4.9% of N is categorized as the HCP, at N_c = 0. Then, 50.2% of N is categorized as the structured coordination at N_c = 0, even when the BAA and VA categorize most of N as Other. ${f}_{{\rm{All}}}\left({{\ell }}_{{\rm{RMSD}}}\right)$ and ${f}_{{\rm{SC}}}\left({{\ell }}_{{\rm{RMSD}}}\right)$ seem to follow Gaussian distributions at N_c = 0, respectively. The deviations of ${f}_{{\rm{All}}}\left({{\ell }}_{{\rm{RMSD}}}\right)$ from the Gaussian become marked, as N_c increases. For example, ${f}_{{\rm{All}}}\left({{\ell }}_{{\rm{RMSD}}}\right)$ has fat tailed distribution at N_c = 10M and has markedly asymmetric form, when 30M≤N_c. In particular, the fat tailed regime of ${f}_{{\rm{All}}}\left({{\ell }}_{{\rm{RMSD}}}\right)$ at ${{\ell }}_{{\rm{RMSD}}}\leqslant 0.3$ is caused by the plateau distribution of ${f}_{{\rm{FCC}}}\left({{\ell }}_{{\rm{RMSD}}}\right)$ at ${{\ell }}_{{\rm{RMSD}}}\leqslant 0.3$ and bimodal distribution of ${f}_{{\rm{BCC}}}\left({{\ell }}_{{\rm{RMSD}}}\right)$ at ${{\ell }}_{{\rm{RMSD}}}\leqslant 0.3$, as shown in figures 44, 45 and 46. Similarly, asymmetry of ${f}_{{\rm{SC}}}\left({{\ell }}_{{\rm{RMSD}}}\right)$ becomes marked, when 10M≤N_c. Other significant result obtained using the PTM is time-evolutions of ${f}_{{\rm{HCP}}}\left({{\ell }}_{{\rm{RMSD}}}\right)$. ${f}_{{\rm{HCP}}}\left({{\ell }}_{{\rm{RMSD}}}\right)$ approaches $\exp \left(-a| {{\ell }}_{{\rm{RMSD}}}-0.27| \right)$ ($a\in {{\mathbb{R}}}_{+}$), as N_c increases. Therefore, most of granular particles, which are categorized as the HCP, deviates from the HCP template with the constant distance 0.27 and such ${f}_{{\rm{HCP}}}\left({{\ell }}_{{\rm{RMSD}}}\right)$ decreases exponentially, as ${{\ell }}_{{\rm{RMSD}}}$ deviates from 0.27.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, the left-top frame of figure 47 shows time-evolution of the fraction of categories XX (:=Other, FCC, HCP, BCC, ICO and SC). The fraction of the HCP increases and becomes comparable with that of the SC, as N_c increases. The temporal change in the fraction of Other is not marked in comparison of that obtained using the BAA and VA. The bottom-left and bottom-right frames show good fits of fractions of categories Other and SC with the inverse-power-law functions, namely, $56{\left({N}_{c}/M+4\right)}^{-0.06}$ and $120{\left({N}_{c}/M+4\right)}^{-0.5}$, respectively. The deviations from the inverse-power-law functions become marked at N_c ≥ 30M in the case of the SC.

Zoom In Zoom Out Reset image size