# Hierarchical structures of amorphous solids characterized by persistent homology

Edited by Giorgio Parisi, University of Rome, Rome, Italy, and approved April 22, 2016 (received for review October 22, 2015)

## Significance

Persistent homology is an emerging mathematical concept for characterizing shapes of data. In particular, it provides a tool called the persistence diagram that extracts multiscale topological features such as rings and cavities embedded in atomic configurations. This article presents a unified method using persistence diagrams for studying the geometry of atomic configurations in amorphous solids. The method highlights hierarchical structures that conventional techniques could not have treated appropriately.

## Abstract

This article proposes a topological method that extracts hierarchical structures of various amorphous solids. The method is based on the persistence diagram (PD), a mathematical tool for capturing shapes of multiscale data. The input to the PDs is given by an atomic configuration and the output is expressed as 2D histograms. Then, specific distributions such as curves and islands in the PDs identify meaningful shape characteristics of the atomic configuration. Although the method can be applied to a wide variety of disordered systems, it is applied here to silica glass, the Lennard-Jones system, and Cu-Zr metallic glass as standard examples of continuous random network and random packing structures. In silica glass, the method classified the atomic rings as short-range and medium-range orders and unveiled hierarchical ring structures among them. These detailed geometric characterizations clarified a real space origin of the first sharp diffraction peak and also indicated that PDs contain information on elastic response. Even in the Lennard-Jones system and Cu-Zr metallic glass, the hierarchical structures in the atomic configurations were derived in a similar way using PDs, although the glass structures and properties substantially differ from silica glass. These results suggest that the PDs provide a unified method that extracts greater depth of geometric information in amorphous solids than conventional methods.

### Sign up for PNAS alerts.

Get alerts for new articles, or get an alert when an article is cited.

The atomic configurations of amorphous solids are difficult to characterize. Because they have no periodicity as found in crystalline solids, only local structures have been analyzed in detail. Although short-range order (SRO) defined by the nearest neighbor is thoroughly studied, it is not sufficient to fully understand the atomic structures of amorphous solids. Therefore, medium-range order (MRO) has been discussed to properly characterize amorphous solids (1–3). Many experimental and simulation studies (4–7) have suggested signatures of MRO such as a first sharp diffraction peak (FSDP) in the structure factor of the continuous random network structure, and a split second peak in the radial distribution function of the random packing structure. However, in contrast to SRO, the geometric interpretation of MRO and the hierarchical structures among different ranges are not yet clear.

Among the available methods, the distributions of bond angle and dihedral angle are often used to identify the geometry beyond the scale of SRO. They cannot, however, provide a complete description of MRO because they only deal with the atomic configuration up to the third nearest neighbors. Alternatively, ring statistics are also applied as a conventional combinatorial topological method (2, 8, 9). However, this method is applicable only for the continuous random network or crystalline structures, and furthermore it cannot describe length scale. Therefore, methodologies that precisely characterize hierarchical structures beyond SRO and are applicable to a wide variety of amorphous solids are highly desired.

In recent years, topological data analysis (10, 11) has rapidly grown and has provided several tools for studying multiscale data arising in physical and biological fields (11–16). A particularly important tool in the topological data analysis is persistence diagram (PD), a visualization of persistent homology as a 2D histogram (e.g., see Fig. 2). The input to the PD is given by an atomic configuration with scale parameters, and the output consists of various multiscale information about topological features such as rings and cavities embedded in the atomic configuration. Here, the atomic configurations are generated by molecular dynamics simulations in this article. Importantly, in contrast to other topological tools, PDs not only count topological features but also provide the scales of these features. Hence, PDs can be used to classify topological features by their scales and clarify geometric relationships among them; this is presumably the most desired function for deeper analysis of amorphous structures.

This article proposes a method using PDs for various amorphous solids in a unified framework. The method is applied to atomic configurations and enables one to study hierarchical geometry embedded in amorphous structures that cannot be treated by conventional methods. We first applied the method to silica glass as an example of the continuous random network structure and obtained the following results. (

*i*) We found three characteristic curves in the PD of silica glass. These curves classify the SRO rings in the ${\text{SiO}}_{4}$ tetrahedra and the MRO rings constructed by those tetrahedra. Furthermore, a hierarchical relationship among the SRO and MRO rings was elucidated. (*ii*) The PD reproduced the wavelength of the FSDP and clarified a real space origin of the FSDP. (*iii*) Each curve in the PD represents a geometric constraint on the ring shapes and, as an example, an MRO constraint on rings consisting of three oxygen atoms was explicitly derived as a surface in a parameter space of the triangles. Moreover, we verified that these curves are preserved under strain, indicating that the PD properly encodes the material property of elastic response. Next, as examples of the random packing structure, the Lennard-Jones (LJ) system and Cu-Zr metallic glass were studied by the PDs, and we clarified the following. (*iv*) These amorphous solids were also characterized well by the distributions of curves and islands in the PDs. (*v*) In the LJ system, the global connectivity of dense packing regions was revealed by dualizing octahedral arrangements. (*vi*) In Cu-Zr alloys, we found that the pair-distribution function defined by the octahedral region in the PD shows the split second peak. Furthermore, a relationship between the hierarchical ring structure and high glass-forming ability was discovered in Cu-Zr alloys.## PDs of Atomic Configuration

The input to PDs is a pair $\mathcal{A}=\left(Q,R\right)$ of an atomic configuration $Q=\left({\mathit{x}}_{1},\dots ,{\mathit{x}}_{N}\right)$ and a parameter set $R=\left({r}_{1},\dots ,{r}_{N}\right)$. Here, ${\mathit{x}}_{i}$ and ${r}_{i}$ are the position in ${\mathrm{\mathbb{R}}}^{3}$ and the input radius for the

*i*th atom, respectively. To characterize the multiscale properties in*Q*, we introduce a parameter*α*, which controls resolution, and generate a family of atomic balls ${B}_{i}\left(\alpha \right)=\left\{\mathit{x}\in {\mathrm{\mathbb{R}}}^{3}|\Vert \mathit{x}-{\mathit{x}}_{i}\Vert \le {r}_{i}\left(\alpha \right)\right\}$ having the radius ${r}_{i}\left(\alpha \right)=\sqrt{\alpha +{r}_{i}^{2}}$. We vary the radii of the atomic balls by*α*and detect rings and cavities at each*α*, where $\alpha \ge {\alpha}_{\text{min}}\u2254-\text{min}\left\{{r}_{1}^{2},\dots ,{r}_{N}^{2}\right\}$.Let ${c}_{k}$ be a ring or cavity constructed in the atomic ball model $B\left(\alpha \right)={\mathrm{\cup}}_{i=1}^{N}{B}_{i}\left(\alpha \right)$ at a parameter

*α*. To be more precise, a ring [respectively (resp.) cavity] here means a generator of the homology ${H}_{1}\left(B\left(\alpha \right)\right)$ [resp. ${H}_{2}\left(B\left(\alpha \right)\right)$] with a field coefficient (17). Then, we observe that there is a value $\alpha ={b}_{k}$ (resp. $\alpha ={d}_{k}$) at which ${c}_{k}$ first appears (resp. disappears) in the atomic ball model. The values ${b}_{k}$ and ${d}_{k}$ are called the birth and death scales of ${c}_{k}$, respectively. The collection of all of the $\left({b}_{k},{d}_{k}\right)\in {\mathrm{\mathbb{R}}}^{2}$ of rings (resp. cavities) is the PD denoted by ${D}_{1}\left(\mathcal{A}\right)$ [resp. ${D}_{2}\left(\mathcal{A}\right)$] for $\mathcal{A}$ (Fig. 1). It follows from the structure theorem of persistent homology (11) that the PDs are uniquely defined from the input. From this construction, $\left({b}_{k},{d}_{k}\right)$ encodes certain scales of each ${c}_{k}$. For example, in ${D}_{1}\left(\mathcal{A}\right)$, ${b}_{k}$ indicates the maximum distance between two adjacent atoms in the ring ${c}_{k}$, whereas ${d}_{k}$ indicates the size of ${c}_{k}$.Fig. 1.

In this article, our basic strategy is that we transform a complicated atomic configuration into PDs and try to identify meaningful shape information from specific distributions such as curves or islands in the PDs. Namely, we reconstruct characteristic atomic subsets from each distribution. To this aim, we compute the optimal cycle for each point $\left({b}_{k},{d}_{k}\right)\in {D}_{\ell}\left(\mathcal{A}\right)$ on the distribution. Mathematically speaking, for a given homology generator ${c}_{k}$ of $\left({b}_{k},{d}_{k}\right)$, the optimal cycle is obtained by solving a minimizing problem in the representatives of ${c}_{k}$ under ${\ell}_{1}$-norm (see refs. 18 and 19). Our method combining PDs with optimal cycles provides a tool to study inverse problems of PDs and effectively works in the geometric analysis of glass structures, as we will see shortly.

For a mathematically rigorous introduction of these concepts see

*Supporting Information*or refs. 10 and 11. In this article, the PDs are computed by CGAL (20) and PHAT (21).## PDs for Continuous Random Network Structure

Fig. 2 shows the PDs ${D}_{1}\left({\mathcal{A}}_{\text{liq}}\right)$, ${D}_{1}\left({\mathcal{A}}_{\text{amo}}\right)$, and ${D}_{1}\left({\mathcal{A}}_{\text{cry}}\right)$ of a liquid ${\mathcal{A}}_{\text{liq}}=\left({Q}_{\text{liq}},{R}_{\text{liq}}\right)$, an amorphous ${\mathcal{A}}_{\text{amo}}=\left({Q}_{\text{amo}},{R}_{\text{amo}}\right)$, and a crystalline ${\mathcal{A}}_{\text{cry}}=\left({Q}_{\text{cry}},{R}_{\text{cry}}\right)$ state of silica, respectively. Here, the horizontal and vertical axes are the birth and death scales, respectively, and the multiplicity of the PDs is plotted on a logarithmic scale. The configurations ${Q}_{\text{liq}}$, ${Q}_{\text{amo}}$, and ${Q}_{\text{cry}}$ are acquired by molecular dynamics simulations using the Beest–Kramer–Santen (BKS) model (22). The input radii

*R*are set to be ${r}_{\text{O}}=1.275$ Å and ${r}_{\text{Si}}=0.375$ Å for each type of the atom (O or Si), which were determined from the first peak positions of the partial radial distribution functions of the amorphous configuration ${Q}_{\text{amo}}$. The details of the molecular dynamics simulations and input radii are given in*Supporting Information*.Fig. 2.

We discovered that the PDs in Fig. 2 distinguish these three states. The liquid, amorphous, and crystalline states are characterized by planar (2-dim), curvilinear (1-dim), and island (0-dim) regions of the distributions, respectively. Here, the 0 and 2 dimensionality of the PDs result from the periodic and random atomic configurations of the crystalline and liquid states, respectively. Furthermore, we emphasize that the presence of the curves in ${D}_{1}\left({\mathcal{A}}_{\text{amo}}\right)$ clearly distinguishes the amorphous state from the others. This implies that specific geometric features of the rings generating the curves in ${D}_{1}\left({\mathcal{A}}_{\text{amo}}\right)$ play a significant role for elucidating amorphous states.

As shown in Fig. 2, ${D}_{1}\left({\mathcal{A}}_{\text{amo}}\right)$ contains three characteristic curves ${C}_{\text{P}}$, ${C}_{\text{T}}$, and ${C}_{\text{O}}$ and one band region ${B}_{\text{O}}$, which are precisely characterized by using the invariance property with respect to the initial radius (15). These particular distributions, especially ${C}_{\text{O}}$, start to become isolated near the glass transition temperature $T={T}_{\text{g}}$ (

*Supporting Information*). Through further analysis of the persistent homology using optimal cycles, we found the following three geometric characterizations. (*i*) The rings on ${C}_{\text{P}}$ generate secondary rings on ${C}_{\text{T}},{C}_{\text{O}}$, and ${B}_{\text{O}}$ (P is named after “primary”). That is, by increasing the parameter*α*, each ring on ${C}_{\text{P}}$ becomes thicker and starts to create new rings by pinching itself, and these newly generated rings appear on ${C}_{\text{T}},{C}_{\text{O}}$, and ${B}_{\text{O}}$ (a schematic illustration of the pinching process is described in Fig. 1). This indicates a hierarchical structure from ${C}_{\text{P}}$ to ${C}_{\text{T}},{C}_{\text{O}}$, and ${B}_{\text{O}}$ in the continuous random network. An example of the rings in this hierarchical relationship is depicted in the inset of ${D}_{1}\left({\mathcal{A}}_{\text{amo}}\right)$. (*ii*) The rings on ${C}_{\text{T}}$ are constructed by tetrahedra consisting of four oxygen atoms at the vertices with one silicon atom at the center. (*iii*) The rings on ${C}_{\text{O}}$ and ${B}_{\text{O}}$ are constructed only by the oxygen atoms (three and more, respectively). Recalling that the death scale indicates the size of rings, the rings on ${C}_{\text{T}}$ are classified as SRO, whereas those on ${C}_{\text{P}}$, ${C}_{\text{O}}$, and ${B}_{\text{O}}$ are classified as MRO.## Decomposition of FSDP

The FSDP observed in the structure factor $S\left(q\right)$ ($q\sim 1.5$–2 ${\mathrm{\AA}}^{-1}$) has been supposed to be a signature of MRO, but its real space origin is still controversial (4). Here, we found that the distributions ${C}_{\text{P}}$, ${C}_{\text{O}}$, and ${B}_{\text{O}}$ of the MRO rings reproduce the

*q*values of the FSDP fairly well. Moreover, we classified the MRO rings as a real space origin of the FSDP.We first note that the death scales of the rings on ${C}_{\text{P}}$, ${C}_{\text{O}}$, and ${B}_{\text{O}}$ are determined only by the oxygen atoms, and this is directly verified by the PD computation. In addition, recall that

*α*is the parameter controlling the radius ${r}_{i}\left(\alpha \right)=\sqrt{\alpha +{r}_{i}^{2}}$ of the*i*th atomic ball ${B}_{i}\left(\alpha \right)$, and the death scale $\alpha ={d}_{k}$ indicates the size of the individual ring ${c}_{k}$. More precisely, the ring ${c}_{k}$ disappears at $\alpha ={d}_{k}$ by being covered up in the atomic ball model ${\cup}_{i=1}^{N}{B}_{i}\left(\alpha \right)$. Hence $\ell \left({d}_{k}\right)=2\sqrt{{d}_{k}+{r}_{\text{O}}^{2}}$ measures the size of ${c}_{k}$ on ${C}_{\text{P}}$, ${C}_{\text{O}}$, and ${B}_{\text{O}}$, where ${r}_{\text{O}}$ is the input radius of the oxygen atom.From the aforementioned argument, we define a distributionwhere ${A}_{1}={C}_{\text{P}}\cup {C}_{\text{O}}\cup {B}_{\text{O}}$, ${A}_{2}={C}_{\text{P}}$, ${A}_{3}={C}_{\text{O}}$, and ${A}_{4}={B}_{\text{O}}$, and $\left|{A}_{i}\right|$ is the number of the elements in ${A}_{i}$. Here,

$${M}_{{A}_{i}}\left(q\right)=\frac{1}{\left|{A}_{i}\right|}{\displaystyle \sum _{\left({b}_{k},{d}_{k}\right)\in {A}_{i}}\delta \left(q-\frac{2\pi}{\ell \left({d}_{k}\right)}\right)},$$

*δ*is the Dirac delta function, which is used to count the contribution of each MRO ring in*q*-space.Fig. 3 shows the plots of ${M}_{{A}_{i}}\left(q\right)$ and the structure factor $S\left(q\right)$ around the FSDP. We found good agreement between the

*q*values of the FSDP and the peak of ${M}_{{A}_{1}}\left(q\right)$. This implies that the MRO rings composed of ${C}_{\text{O}}$, ${C}_{\text{P}}$, and ${B}_{\text{O}}$ are the real space origin of the FSDP. We also note that the distributions ${M}_{{A}_{3}}\left(q\right),{M}_{{A}_{2}}\left(q\right),{M}_{{A}_{4}}\left(q\right)$ are located on the large, medium, and small*q*values and, hence, the rings in ${C}_{\text{O}}$, ${C}_{\text{P}}$, and ${B}_{\text{O}}$ provide a decomposition of the FSDP into those*q*values, respectively. It should be emphasized that the ${M}_{{A}_{i}}\left(q\right)$ are derived from the configuration of the oxygen atoms only. This shows that the configuration of oxygen atoms plays a significant role in the FSDP. Furthermore, the invariant property (15) of $\ell \left({d}_{k}\right)$ for the rings on ${C}_{\text{P}}$, ${C}_{\text{O}}$, and ${B}_{\text{O}}$ induces that of ${M}_{{A}_{i}}\left(q\right)$ under the choice of the input radius ${r}_{\text{O}}$. This means that our analysis using PDs does not contain any artificial ambiguity of the input radii.Fig. 3.

## Curves and Shape Constraints

The presence of curves ${C}_{\text{P}}$, ${C}_{\text{T}}$, and ${C}_{\text{O}}$ in ${D}_{1}\left({\mathcal{A}}_{\text{amo}}\right)$ clearly distinguishes the amorphous state from the crystalline and liquid states. We emphasize here that this characteristic property shows the constraints on the shapes of the rings induced by the normal directions of these curves.

For example, the shape of a ring on ${C}_{\text{O}}$, which consists of three oxygen atoms (Fig. 4,

*Right*), is determined by specifying the first and second minimum edge lengths ${d}_{1}$ and ${d}_{2}$ (${d}_{1}<{d}_{2}$) and the angle*θ*between them, and hence is realized in a 3D parameter space. Then, the constraint for ${C}_{\text{O}}$ to be the curve requires that these three variables $\left({d}_{1},{d}_{2},\theta \right)$ satisfy a certain relation $f\left({d}_{1},{d}_{2},\theta \right)=0$, and hence provides a restriction on the shape of the O-O-O triangles. Fig. 4 shows a plot of $\left({d}_{1},{d}_{2},\theta \right)$ for the O-O-O triangles on ${C}_{\text{O}}$ in ${D}_{1}\left({\mathcal{A}}_{\text{amo}}\right)$, and we found a surface corresponding to this constraint $f\left({d}_{1},{d}_{2},\theta \right)=0$ in the parameter space. This demonstrates one of the medium-range geometric structures in the amorphous state. It is worth noting that this previously unidentified characterization of MRO of O-O-O triangles in 3D parameter space cannot be derived by separately analyzing the conventional distributions of lengths or angles, because each of them can only deal with the single parameter.Fig. 4.

## Response Under Strain

The presence of the curves also indicates variations in the shapes of the rings. That is, by following each curve along its tangential direction and studying its rings, we can observe the deformation of the rings. It is reasonable to suppose that these variations are due to thermal fluctuations, and hence the deformations of the ring configurations along the curves are probably softer than those in the normal directions. Consequently, the response under strain is expected to follow the same shape constraints.

To verify this mechanical response of PDs, we performed simulations of isotropic compression for our amorphous state and computed the PD of the strained state. The strain is set to be $1\mathrm{\%}$ of the original volume, which is sufficiently small to satisfy a linear response relation with the stress. The three panels on the top in Fig. 5 show the contours of the histograms restricted on ${C}_{\text{P}}$, ${C}_{\text{O}}$, and ${C}_{\text{T}}$ for the original configuration ${Q}_{\text{amo}}$ of the amorphous state (black), for the strained state ${Q}_{\text{amo}}^{\text{strain}}$ (red), and for the artificial $1\mathrm{\%}$ linear shrink ${Q}_{\text{amo}}^{\text{linear}}$ of the original coordinates (blue), respectively. The contours of ${C}_{\text{O}}$ are depicted on the coordinates along the tangential direction. The bottom three panels in Fig. 5 show the projections of the histograms on their normal axes.

Fig. 5.

Fig. 5 shows that the contours of the strained state shift along the original curves ${C}_{\text{P}}$, ${C}_{\text{O}}$, and ${C}_{\text{T}}$ in ${D}_{1}\left({\mathcal{A}}_{\text{amo}}\right)$ (i.e., downward in ${C}_{\text{P}}$, leftward in ${C}_{\text{O}}$, and almost fixed in ${C}_{\text{T}}$, respectively). This is in contrast to a linear shrink, in which the contours simply move in the direction of decreasing both birth and death scales because of the uniform shrink of the system size. This strongly suggests that the configuration ${Q}_{\text{amo}}^{\text{strain}}$ reflects the shape constraints and supports the expected mechanical response of PDs. We also note from the figures of ${Q}_{\text{amo}}^{\text{strain}}$ that the contour of ${C}_{\text{T}}$ is almost fixed compared with those of ${C}_{\text{P}}$ and ${C}_{\text{O}}$. Because ${C}_{\text{P}}$ and ${C}_{\text{O}}$ (${C}_{\text{T}}$) are classified as MRO (SRO), this implies that MRO is mechanically softer than SRO. Here, we remark that we can observe a similar behavior of the PDs for a shear deformation.

We have revealed that PDs encode information about elastic response, similar to how the radial distribution function encodes volume compressibility (23). It should also be emphasized that the curves or islands appear in the PDs of only the solids. This evidence suggests that these isolated distributions are related to the rigidity of the materials. This hypothesis follows from the fact that isolations represent geometric constrains reflecting mechanical responses. Future numerical and theoretical studies to unravel this relationship would be of great value.

## PDs for Random Packing Structures

We next study the geometry of amorphous states close to random packing structures. In this case, we found that both ${D}_{1}$ and ${D}_{2}$ capture characteristic features of the amorphous structures. Fig. 6 shows the PDs ${D}_{i}\left({\mathcal{A}}_{\text{cry}}^{\text{LJ}}\right)$ and ${D}_{i}\left({\mathcal{A}}_{\text{amo}}^{\text{LJ}}\right)$ for $i=1,2$ of the monodisperse LJ system in the crystalline and the amorphous states, respectively. The input radii $r={r}_{1}=\dots ={r}_{N}$ in

*R*are set to be zero, because changing*r*only causes translations of the PDs for the single component system. The details of the simulation are explained in*Supporting Information*.Fig. 6.

Similar to the case of silica, the crystalline structure is characterized by the island distributions in the PDs (top panels in Fig. 6). They correspond to the regular triangles in ${D}_{1}\left({\mathcal{A}}_{\text{cry}}^{\text{LJ}}\right)$ and the regular octahedra, the regular tetrahedra, and the quartoctahedra (24) in ${D}_{2}\left({\mathcal{A}}_{\text{cry}}^{\text{LJ}}\right)$ of the face-centered cubic (FCC) configuration. For the amorphous structure, the curves in ${D}_{1}\left({\mathcal{A}}_{\text{amo}}^{\text{LJ}}\right)$ and ${D}_{2}\left({\mathcal{A}}_{\text{amo}}^{\text{LJ}}\right)$ represent variations of triangles and tetrahedra, respectively. We also note that the isolation of the octahedral distribution is preserved well even in ${D}_{2}\left({\mathcal{A}}_{\text{amo}}^{\text{LJ}}\right)$. Its peak is separated from the curve of the tetrahedra (bottom right panel in Fig. 6), demonstrating the quantitative classification into two typical local structures of different-sized cavities.

In random packings, it is known that the atomic configuration can be divided into dense packing regions built from tetrahedra and the complement that patches those regions together (25). In particular, the network structure of the dense packing regions characterizes the global connectivity beyond MRO, which has not yet been investigated in detail. Note that ${D}_{2}\left({\mathcal{A}}_{\text{amo}}^{\text{LJ}}\right)$ clearly separates the dense packing regions as the deformation curve of the tetrahedra and identifies the octahedral island as the main component of the complement structure. From the Alexander duality in ${\mathrm{\mathbb{R}}}^{3}$ (e.g., ref. 17), the connectivity of the dense packing regions can be studied by the cavities of the complement. Namely, we first extract all octahedra from the octahedral island in ${D}_{2}\left({\mathcal{A}}_{\text{amo}}^{\text{LJ}}\right)$ and construct a new set of points ${\mathcal{A}}_{\text{oct}}^{\text{LJ}}=\left({Q}_{\text{oct}}^{\text{LJ}},{R}_{\text{oct}}^{\text{LJ}}\right)$ by putting a point at the center of each octahedron (Fig. 7,

*Left*). Here, we set ${R}_{\text{oct}}^{\text{LJ}}$ to be zero, because the SD of the octahedral sizes is very small.Fig. 7.

Fig. 7,

*Right*shows ${D}_{2}\left({\mathcal{A}}_{\text{oct}}^{\text{LJ}}\right)$, and we clearly observe that almost all cavities are located close to the diagonal. This means that the octahedral distribution rarely generates persistent cavities, and hence by the duality we can conclude that the dense packing regions mostly construct a giant connected network. It should be remarked that the dual treatment is much easier and computationally more efficient than directly studying the intricate huge network structures. We also emphasize that the analysis here of treating the octahedra as a new input is an iterative use of the PD method and can be an effective approach for studying multiscale geometry.As an example of multicomponent systems, we also studied metallic glasses composed of Cu and Zr (26), in particular, focusing on ${\text{Cu}}_{50}{\text{Zr}}_{50}$ and ${\text{Cu}}_{15}{\text{Zr}}_{85}$, which display the different glass-forming ability. The PDs for these alloys are shown in Fig. 8. The input radii for Cu and Zr are set to be 1.30 Å and 1.55 Å, respectively, for the multicomponent PDs (Fig. 8,

*Top*). These values are obtained by the same procedure as for the silica. Even in the multicomponent system, the PDs basically show similar behaviors to those of the LJ system. Specifically, the island distribution corresponding to octahedra appears in ${D}_{2}\left({\mathcal{A}}_{\text{amo}}^{{\text{Cu}}_{50}{\text{Zr}}_{50}}\right)$ and the characteristic curves are also observed in ${D}_{1}\left({\mathcal{A}}_{\text{amo}}^{{\text{Cu}}_{50}{\text{Zr}}_{50}}\right)$ and ${D}_{2}\left({\mathcal{A}}_{\text{amo}}^{{\text{Cu}}_{50}{\text{Zr}}_{50}}\right)$.Fig. 8.

In the random packing structure, the split second peak of the radial distribution function has been supposed to be a signature of MRO (27). The shaded region of the radial distribution function in Fig. 9,

*Bottom*shows the split second peak of ${\text{Cu}}_{50}{\text{Zr}}_{50}$. Then, we found that the pair-distance distribution of the atoms of generators in $B\u2254\left[0.92,1.05\right]\times \left[1.76,2.01\right]\subset {D}_{2}\left({\mathcal{A}}_{\text{amo}}^{{\text{Cu}}_{50}{\text{Zr}}_{50}}\right)$ also shows a clear splitting of the distribution (black line in Fig. 9,*Top*) in the same length scale. Here,*B*is chosen to be a region around the octahedral distribution. Meanwhile the pair-distance distribution of generators other than*B*shows a slight change there (pink line in Fig. 9). This means that the generators around the octahedral distribution play a significant role for the split second peak. Therefore, similar to the FSDP in the silica, this result demonstrates that the PDs classify the length scale of MRO from other scales.Fig. 9.

We also studied the PDs of only the Zr component. The PD ${D}_{1}$ for Zr in ${\text{Cu}}_{50}{\text{Zr}}_{50}$ (Fig. 8,

*Middle Left*) represents the existence of a hierarchical MRO structure similar to that of the silica in Fig. 2, whereas the PD ${D}_{1}$ for Zr in ${\text{Cu}}_{15}{\text{Zr}}_{85}$ (Fig. 8,*Bottom Left*) does not show any hierarchical curves. An example of the hierarchical rings in ${\text{Cu}}_{50}{\text{Zr}}_{50}$ is depicted as the insets in the PD.We here remark that ${\text{Cu}}_{50}{\text{Zr}}_{50}$ has higher glass-forming ability than ${\text{Cu}}_{15}{\text{Zr}}_{85}$ (28). Relationships between the glass-forming ability and the geometric structure are now actively studied (e.g., ref. 29). Then, this result suggests another possibility that the presence of the hierarchical MRO structure extracted from PDs is also related to the glass-forming ability of the alloy. In view of the results that the hierarchical MRO structure in PDs plays an important role for characterizing the glass states in ${\text{SiO}}_{2}$, this statement seems to be reasonable. To understand the geometric and topological formation of the glass, it will be a great challenge to clarify this new perspective.

## Conclusion

We have presented that PD is a powerful tool for geometric characterizations of various amorphous solids in the short, medium, and even further ranges. In this work, we have addressed two different types of amorphous systems: continuous random network and random packing structures. Both types of amorphous systems are characterized well by the existence of the curve and island distributions in the PDs. These specific distributions characterize the shapes of rings and cavities in multiranges, and the analysis using optimal cycles explicitly captures hierarchical structures of these shapes. We have shown that these shape characteristics successfully reproduce the FSDP for the continuous random network and the split second peak for the random packing and provide further geometric insights to them. Furthermore, the global connectivity of dense packing regions in the LJ system is revealed by the iterative application of the PD. For the binary random packing of the Cu-Zr metallic glass, we have also shown that the presence of the hierarchical MRO rings in the single component suggests the relationship with the glass-forming ability.

The methodology presented here can be applied to a wide variety of disordered systems and enables one to survey the geometric features and constraints in seemingly random configurations. Furthermore, because we investigated the mechanical response of the PDs, studying dynamical properties of materials using the PD method would be of great importance to understand the relationship between hierarchical structures and mechanical properties. We believe that further developments and applications of topological data analysis will accelerate the understanding of amorphous solids.

## Computational Homology

Here we summarize the computational homological tools used in this article. For further mathematical and computational details, please refer to ref. 10.

Given a pair $\mathcal{A}=\left(Q,R\right)$ of an atomic configuration $Q=\left({\mathit{x}}_{1},{\mathit{x}}_{2}\dots ,{\mathit{x}}_{N}\right)$ and a set $R=\left({r}_{1},{r}_{2},\dots ,{r}_{N}\right)$ of input radii, the persistence diagram captures hierarchical structures of the topological features such as rings and cavities in a one-parameter family of alpha shapes $A\left(Q,R;\alpha \right)$. Here, the alpha shape $A\left(Q,R;\alpha \right)$ is constructed as follows. Let ${\mathrm{\mathbb{R}}}^{3}={\cup}_{i=1}^{N}{V}_{i}$ be the Voronoi decomposition for

*Q*, where ${V}_{i}$ is the Voronoi cell of ${\mathit{x}}_{i}$. For each parameter*α*, we assign a ball ${B}_{i}\left(\alpha \right)=\left\{\mathit{x}\in {\mathrm{\mathbb{R}}}^{3}|\Vert \mathit{x}-{\mathit{x}}_{i}\Vert \le {r}_{i}\left(\alpha \right)\right\}$ with the radius ${r}_{i}\left(\alpha \right)=\sqrt{\alpha +{r}_{i}^{2}}$ for the*i*th atom. Then, the alpha shape $A\left(Q,R;\alpha \right)$ is defined by the dual of the cut balls ${C}_{i}\left(\alpha \right)={V}_{i}\mathrm{\cap}{B}_{i}\left(\alpha \right)$, $i=1,\dots ,N$. Namely, this is a polyhedron with the vertex set*Q*such that an edge $\left|{\mathit{x}}_{i}{\mathit{x}}_{j}\right|$ (triangle and tetrahedron, respectively) is assigned to $A\left(Q,R;\alpha \right)$ if and only if the intersection ${C}_{i}\left(\alpha \right)\mathrm{\cap}{C}_{j}\left(\alpha \right)$ of the corresponding double (triple and quadruple, respectively) cut balls is nonempty (Fig. S1). In this article, the alpha shapes were computed using CGAL (20).Fig. S1.

One of the important properties of the alpha shape is that the atomic ball model $B\left(\alpha \right)={\mathrm{\cup}}_{i=1}^{N}{B}_{i}\left(\alpha \right)$ and the alpha shape $A\left(Q,R;\alpha \right)$ can be continuously deformed to each other. Hence, there is no loss of topological information in the use of the alpha shape model, which is much easier to analyze computationally than the atomic ball model. Furthermore, we can trace the hierarchical structures in $B\left(\alpha \right)$ by changing the parameter $\alpha \ge {\alpha}_{\text{min}}$, where ${\alpha}_{\text{min}}=-\text{min}\left\{{r}_{1}^{2},\dots {r}_{N}^{2}\right\}$. In other words, the parameter

*α*controls the resolution from fine ($\alpha \approx {\alpha}_{\text{min}}$) to crude ($\alpha \gg {\alpha}_{\text{min}}$) geometric features.Fig. S2,

*Left*shows a schematic illustration of atomic ball models $B\left(\alpha \right)$ with the rings on ${C}_{\text{P}}$ (red), ${C}_{\text{T}}$ (blue), ${C}_{\text{O}}$ (yellow), and ${B}_{\text{O}}$ (green). Here, the large (small) balls correspond to oxygen (silicon). At each ${\alpha}_{i}$ ($i=2,3,4,5$), a new ring appears, and the dashed rings express the rings that have disappeared. The hierarchical ring structures are observed at ${\alpha}_{3}$, ${\alpha}_{4}$, and ${\alpha}_{5}$, at which new rings are generated by pinching the primary ring.Fig. S2.

The persistence diagram is a 2D histogram recording the birth and death scales of the topological features in a one-parameter family ${\left\{B\left(\alpha \right)\right\}}_{\alpha}$ of atomic ball models. For example, let

*c*be the yellow ring in Fig. S2. We observe that*c*first appears at $\alpha ={\alpha}_{4}$ and first disappears at some value ${\alpha}_{\text{*}}$, where ${\alpha}_{4}<{\alpha}_{\text{*}}<{\alpha}_{5}$. Then, the birth and death scales*b*and*d*, respectively, of the ring*c*are determined by $b={\alpha}_{4}$ and $d={\alpha}_{*}$. The histogram of the birth and death pairs $\left(b,d\right)\in {\mathrm{\mathbb{R}}}^{2}$ for all of the rings appearing in the one-parameter family ${\left\{B\left(\alpha \right)\right\}}_{\alpha}$ is called the persistence diagram ${D}_{1}\left(\mathcal{A}\right)$. Fig. S2,*Right*shows the persistence diagram of the left. From the definition of the persistence diagram, the birth scale indicates the maximum distance between adjacent atoms in the ring, whereas the death scale indicates the size of the ring. It should also be mentioned that persistence diagrams ${D}_{i}\left(\mathcal{A}\right)$ can be similarly defined for any nonnegative integer $i\in \left\{0,1,2,\dots \right\}$, and they capture higher-dimensional topological features (e.g., cavities for $i=2$). The persistence diagrams in this article were computed by PHAT (21).## MD Simulation for the Silica System

The atomic configurations of silica in a liquid ${Q}_{\text{liq}}$, an amorphous solid ${Q}_{\text{amo}}$, and a crystalline solid ${Q}_{\text{cry}}$ state were generated by the MD simulation. The system was composed of 2,700 silicon atoms and 5,400 oxygen atoms interacting via the BKS potential (22), and the 24-6 LJ potential is assigned to correct for the unphysical behavior during equilibration in the liquid state (30). The masses of the silicon and oxygen atoms were 27.98 and 15.99 g/mol, respectively. Despite its simplicity, the BKS model reasonably reproduces the static and dynamic properties of amorphous silica (31).

Starting from the beta-cristobalite structure, the crystalline configuration ${Q}_{\text{cry}}$ was obtained after the equilibration for 5 ps at 1 bar and 10 K. The MD simulation was performed with Parrinello–Raman barostat and Nosé–Hoover thermostat with the damping parameters ${\tau}_{\text{bar}}=0.5$ ps and ${\tau}_{\text{th}}=0.8$ ps, respectively. The time step was set to be 1.0 fs.

The liquid configuration ${Q}_{\text{liq}}$ was obtained after the equilibrating the system from a random initial configuration for 50 ps at 1 bar and 7,000 K. Here, the MD simulation was performed with the Andersen barostat and the Nosé–Hoover thermostat.

The amorphous configuration ${Q}_{\text{amo}}$ was obtained after the cooling MD simulation (31). Starting from the liquid configuration at the temperature 7,000 K, the heat-bath temperature was decreased linearly to the final temperature 10 K with $-dT/dt={10}^{12}$ K/s. Throughout the cooling process the pressure was kept at 1 bar and we used the Nosé–Hoover thermostat with ${\tau}_{\text{th}}=1.0$ ps. In the high-temperatures region 1,500 K $\le T\le $ 6,000 K, we used the Anderson barostat with ${\tau}_{\text{bar}}=1.0$ ps. In the low-temperatures region 10 K $\le T\le $ 1,500 K, we used the Parrinello–Raman barostat with ${\tau}_{\text{bar}}=0.5$ ps. Then, ${Q}_{\text{amo}}$ was obtained as the final configuration of the cooling process. We tested five cooling rates, $-dT/dt={10}^{15},{10}^{14},{10}^{13},{10}^{12}$, and ${10}^{11}$ K/s, and checked that these choices do not affect our results qualitatively.

The set $R=\left({r}_{1},\cdots ,{r}_{N}\right)$ of input radii can be adjusted for each

*i*th atom ($i=1,\dots ,N$). In our analysis, we chose the radii ${r}_{\text{Si}}$ and ${r}_{\text{O}}$ of silicon and oxygen using the first peak positions of the partial radial distribution functions of ${Q}_{\text{amo}}$. The first peaks for SiO, OO, and SiSi appeared at ${r}_{\text{SiO}}=1.65$ Å, ${r}_{\text{OO}}=2.55$ Å, and ${r}_{\text{SiSi}}=3.15$ Å, respectively. Then, we selected the smallest two peaks and determined ${r}_{\text{Si}}$ and ${r}_{\text{O}}$ by solving ${r}_{\text{Si}}+{r}_{\text{O}}={r}_{\text{SiO}}$ and $2{r}_{\text{O}}={r}_{\text{OO}}$. This led to ${r}_{\text{O}}=1.275$ Å and ${r}_{\text{Si}}=0.375$ Å.The glass transition temperature ${T}_{\text{g}}$ and the melting temperature ${T}_{\text{m}}$ are estimated through the enthalpy temperature curves. Starting from beta-cristobalite configuration, the system is heated with a heating rate $dT/dt=0.1$ K/ps from $T=1$ K to $T=\mathrm{7,000}$ K (black curve in Fig. S3). Then, ${T}_{\text{m}}$ is evaluated to be around $\mathrm{3,900}$ K (32). Subsequently, the cooling simulation is performed with a cooling rate $-dT/dt=10.0$ K/ps (red curve). Then the glass transition temperature is evaluated as an intersection point between tangential lines of low (blue line) and high (green line) temperature regions of the curve and found to be ${T}_{\text{g}}\approx \mathrm{2,700}$ K (31). Fig. S4 shows PDs for $\mathrm{2,500}$ K (

*Left*) and $T=\mathrm{3,500}$ K (*Right*), and the curve ${C}_{\text{O}}$ starts to appear below ${T}_{\text{g}}$.Fig. S3.

Fig. S4.

## MD Simulation for the LJ System and the Cu-Zr Alloy

The LJ system is composed of 4,000 monatomic particles interacting via LJ potential $u\left(r\right)=4\epsilon \left({\left(\sigma /r\right)}^{12}-{\left(\sigma /r\right)}^{6}\right)$ with

*ε*,*σ*, and mass equal to unity. The time step is set to be 0.005 in LJ units.Starting from the FCC configuration, ${Q}_{\text{cry}}^{\text{LJ}}$ is obtained after equilibration at temperature $T=0.1$ and pressure $p=1$ with ${\tau}_{\text{th}}=1.0$ and ${\tau}_{\text{bar}}=1.0$. To obtain an amorphous structure ${Q}_{\text{amo}}^{\text{LJ}}$, the equilibration at $T=2$ in the normal liquid phase has been performed for 2,000 steps at the constant number density 1.015 (33) with ${\tau}_{\text{th}}=0.5$. Then, the cooling MD simulations have been performed for several cooling ratios from $-dT/dt={10}^{-5}$ to ${10}^{-1}$. In Fig. S5,

*Left*the enthalpy normalized by the number of particles is described as a function of*T*. As is observed, there is no glass transition in the LJ system. Therefore, we use an inherent structure for the liquid state as ${Q}_{\text{amo}}^{\text{LJ}}$ (black points in Fig. S5,*Left*).Fig. S5.

For the Cu-Zr alloy system, MD simulations have been performed using the embedded-atom method potential (34). The masses of Cu and Zr atoms were set to be 63.54 and 91.22 g/mol, respectively. The number of particles was 16,000, in which there were 8,000 of both copper and zirconium atoms for ${\text{Cu}}_{50}{\text{Zr}}_{50}$, 2,400 were copper atoms, and 13,600 were zirconium atoms for ${\text{Cu}}_{15}{\text{Zr}}_{85}$.

In this case, we can observe the glass transition (Fig. S5,

*Right*) during a cooling simulation with the quench rate ${10}^{13}$ K/s (26) down to $T=1$ K keeping the pressure $p=0$ atm after the equilibration at $T=\mathrm{17,000}$ K. Then, we use the final configuration at $T=1$ K as the amorphous configuration.All MD simulations in this article were performed using LAMMPS (35).

## Acknowledgments

We thank Mingwei Chen, Hajime Tanaka, Masakazu Matsumoto, and Daniel Miles Packwood for valuable discussions and comments. This work was sponsored by World Premier International Research Center Initiative, Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. This work was also supported by Japan Science and Technology Agency (JST) CREST Mathematics Grant 15656429 (to Y.H.), Structural Materials for Innovation Strategic Innovation Promotion Program D72 (Y.H., A.H., and Y.N.), MEXT Coop with Math Program (K.M.), Japan Society for the Promotion of Science (JSPS) Grant 26310205 (to Y.N.), and JSPS Grant 15K13530 and JST PRESTO (to T.N.).

## Supporting Information

Supporting Information (PDF)

Supporting Information

- Download
- 857.47 KB

## References

1

JD Bernal, The Bakerian Lecture, 1962 The structure of liquids.

*Proc R Soc Lond A Math Phys Sci***280**, 299–322 (1964).2

DE Polk, Structural model for amorphous silicon and germanium.

*J Non-Cryst Solids***5**, 365–376 (1971).3

JL Finney, Random packings and the structure of simple liquids. I. The geometry of random close packing.

*Proc R Soc Lond A Math Phys Sci***319**, 479–493 (1970).4

SR Elliott, Origin of the first sharp diffraction peak in the structure factor of covalent glasses.

*Phys Rev Lett***67**, 711–714 (1991).5

SR Elliott

*Physics of Amorphous Materials*(Longman, 2nd Ed, New York, 1990).6

S Susman, et al., Intermediate-range order in permanently densified vitreous SiO2: A neutron-diffraction and molecular-dynamics study.

*Phys Rev B Condens Matter***43**, 1194–1197 (1991).7

GN Greaves, S Sen, Inorganic glasses, glass-forming liquids and amorphizing solids.

*Adv Phys***56**, 1–166 (2007).8

SL Roux, V Petkov, ISAACS – Interactive structure analysis of amorphous and crystalline systems.

*J Appl Cryst***43**, 181–185 (2010).9

WH Zachariasen, The atomic arrangement in glass.

*J Am Chem Soc***54**, 3841–3851 (1932).10

H Edelsbrunner, J Harer,

*Computational Topology: An Introduction*(Am Mathematical Soc, Providence, RI). (2010).11

G Carlsson, Topology and data.

*Bull Am Math Soc***46**, 255–308 (2009).12

M Nicolau, AJ Levine, G Carlsson, Topology based data analysis identifies a subgroup of breast cancers with a unique mutational profile and excellent survival.

*Proc Natl Acad Sci USA***108**, 7265–7270 (2011).13

A Hirata, et al., Geometric frustration of icosahedron in metallic glasses.

*Science***341**, 376–379 (2013).14

H Kurtuldu, K Mischaikow, MF Schatz, Extensive scaling from computational homology and Karhunen-Loève decomposition analysis of Rayleigh-Bénard convection experiments.

*Phys Rev Lett***107**, 034503 (2011).15

T Nakamura, Y Hiraoka, A Hirata, EG Escolar, Y Nishiura, Persistent homology and many-body atomic structure for medium-range order in the glass.

*Nanotechnology***26**, 304001 (2015).16

L Kondic, et al., Topology of force networks in compressed granular media.

*Europhys Lett***97**, 54001 (2012).17

A Hatcher

*Algebraic Topology*(Cambridge Univ Press, Cambridge, UK, 2001).18

T Dey, et al., Optimal homologous cycles, total unimodularity, and linear programming.

*SIAM J Comput***40**, 1026–1044 (2011).19

EG Escolar, Y Hiraoka, Optimal cycles for persistent homology via linear programming.

*Optimization in the Real World –Towards Solving Real-World Optimization Problems*. Mathematics for Industry, eds Fujisawa K, Shinano Y, Waki H (Springer, Tokyo), pp 79–96. (2015).20

CGAL (2016) The Computational Geometry Algorithms Library. Available at www.cgal.org/.

21

PHAT (2014) Persistent Homology Algorithm Toolbox. Available at https://bitbucket.org/phat-code/phat.

22

GJ Kramer, Force fields for silicas and aluminophosphates based on ab initio calculations.

*Phys Rev Lett*; van Beest BW; van Santen RA**64**, 1955–1958 (1990).23

JP Hansen, IR McDonald

*Theory of Simple Liquids*(Academic, 3rd Ed, New York, 2006).24

AV Anikeenko, et al., A novel Delaunay simplex technique for detection of crystalline neuclei in dense packings of spheres.

*Computational Science and Its Applications – ICCSA 2005*. Lecture Notes in Computer Science, eds Gervasi O, Gavrilova ML, Kumar V, Laganà A, Lee HP, Mun Y, Taniar D, Tan CJK (Springer, New York), Vol 3480, pp 816–826. (2005).25

R Yamamoto, M Doyama, The polyhedron and cavity analyses of a structural model of amorphous iron.

*J Phys F Metal Phys***9**, 617–627 (1979).26

SP Pan, et al., Origin of splitting of the second peak in the pair-distribution function for metallic glasses.

*Phys Rev B***84**, 092201 (2011).27

JA Barker, et al., Relaxation of the Bernal model.

*Nature***257**, 120–122 (1975).28

E Ma, Tuning order in disorder.

*Nat Mater***14**, 547–552 (2015).29

YQ Cheng, E Ma, Atomic-level structure and structure-property relationship in metallic glasses.

*Prog Mater Sci***56**, 379–473 (2011).30

Y Guissani, B Guillot, A numerical investigation of the liquid-vapor coexistence curve of silica.

*J Chem Phys***104**, 7633 (1996).31

K Vollmayr, W Kob, K Binder, Cooling-rate effects in amorphous silica: A computer-simulation study.

*Phys Rev B Condens Matter***54**, 15808–15827 (1996).32

TF Soules, et al., Silica molecular dynamic force fields – a practical assessment.

*J Non-Cryst Solids***357**, 1564–1573 (2011).33

G Monaco, S Mossa, Anomalous properties of the acoustic excitations in glasses on the mesoscopic length scale.

*Proc Natl Acad Sci USA***106**, 16907–16912 (2009).34

MI Mendelev, et al., Using atomistic computer simulations to analyze x-ray diffraction data from metallic glasses.

*J Appl Phys***102**, 043501 (2007).35

S Plimpton, Fast parallel algorithms for short-range molecular dynamics.

*J Comput Phys***117**, 1–19 (1995).## Information & Authors

### Information

#### Published in

#### Classifications

#### Copyright

Freely available online through the PNAS open access option.

#### Submission history

**Published online**: June 13, 2016

**Published in issue**: June 28, 2016

#### Keywords

#### Acknowledgments

We thank Mingwei Chen, Hajime Tanaka, Masakazu Matsumoto, and Daniel Miles Packwood for valuable discussions and comments. This work was sponsored by World Premier International Research Center Initiative, Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. This work was also supported by Japan Science and Technology Agency (JST) CREST Mathematics Grant 15656429 (to Y.H.), Structural Materials for Innovation Strategic Innovation Promotion Program D72 (Y.H., A.H., and Y.N.), MEXT Coop with Math Program (K.M.), Japan Society for the Promotion of Science (JSPS) Grant 26310205 (to Y.N.), and JSPS Grant 15K13530 and JST PRESTO (to T.N.).

#### Notes

This article is a PNAS Direct Submission.

### Authors

#### Competing Interests

The authors declare no conflict of interest.

## Metrics & Citations

### Metrics

#### Citation statements

#### Altmetrics

### Citations

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

#### Cited by

Loading...

## View Options

### View options

#### PDF format

Download this article as a PDF file

DOWNLOAD PDF### Get Access

#### Login options

Check if you have access through your login credentials or your institution to get full access on this article.

Personal login Institutional Login#### Recommend to a librarian

Recommend PNAS to a Librarian#### Purchase options

Purchase this article to get full access to it.