Medium-range order in amorphous ices revealed by persistent homology

Despite the amorphous nature of glassy water, x-ray or neutron scattering experiments reveal sharp peaks in the structure factor, indicating the existence of medium-range order (MRO) in the system. However the real space origin of the peaks has yet to be disclosed. Herein, we use a combined approach of molecular dynamics simulations and persistent homology (PH) to investigate two types of glassy water, low-density amorphous (LDA) and high-density amorphous (HDA) ices. We present prominent MRO ring structures in each type of the ices, distinguished by their size and shape as well as the number of their components: MRO rings in HDA are observed smaller, less planar and more membered, compared to those in LDA. The PH-extracted MRO rings successfully reproduce the quantitative features, including the position and width, of the first sharp diffraction peaks in the structure factor, hence suitably serving as the origin of experimental MRO signatures in the amorphous ices. Our study supports that PH is an effective tool to identify hidden MRO in amorphous configurations.


Introduction
Water is likely to crystallize below 0 • C to create ordered structures, which can be of at least 17 different types [1]. In contrast to the polymorphism, solid water can also exist in several distinct metastable disordered states, depending on temperature and pressure. Found in most abundance in the interstellar media [2,3], amorphous ices (or glassy water) were first discovered in an experiment in 1935 [4]. Since then, several experimental routes have been introduced to obtain amorphous ices below the glass transition temperature (120-160 K), such as depositing water vapor on a cold substrate, quenching liquid water or compressing the hexagonal ice all at a sufficiently rapid rate [1,[4][5][6][7]. Among different types of glassy water, two of the most interesting are low-density amorphous (LDA) and high-density amorphous (HDA) ices [8][9][10].
In the long history of our tremendous interest in water, many attempts have been made to study the structure of water and relate it to various properties of water [6,7,[11][12][13][14][15][16][17][18][19][20]. However, it is well-known that in every phase of water the local tetrahedral network (inset figure in figure 1(a)) framed by oxygen atoms and hydrogen bonds (H-bonds) plays an important role in constituting the basic building block of water structure and is regarded to be responsible for water's anomalous behaviors. For instance, the density maximum anomaly of liquid water can be explained by quantifying the degree of tetrahedrality [15], and the wide supercooling region of water can be ascribed to pentagonal rings made of oxygen atoms [16]. The structure of amorphous ices has also been understood in the context of S Supplementary material for this article is available online (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 1 Author to whom any correspondence should be addressed. tetrahedrality. The LDA phase is regarded to maintain locally favored tetrahedral structures, similar to liquid water [17]. The HDA phase, on the other hand, has been understood to be a collapsed and poorly crystalline phase due to the penetration of another water molecule in the second coordination shell toward the inner first shell under pressure [21].
As such, most previous structural analyses of water have been limited to the short-range order (SRO) only. In amorphous water systems, SRO mostly refers to the size comparable to that of the local tetrahedral network frame, which is less than 3.0 Å while medium-range order (MRO) corresponds to the size range beyond that. In contrast to the SRO, our understanding of MRO in amorphous water systems still remains poor today. It has long been observed in x-ray or neutron diffraction experiments of water that the peak intensity in the structure factor increases as the temperature gets lowered, which directly implies an increase of MRO along the glass transition [14,22]. In particular, for a glassy state of water below 120 K, the first sharp diffraction peak (FSDP) in the structure factor supports the existence of strong MRO features embedded in the amorphous system. To the best of our knowledge, however, the structural real space origin of those sharp peaks in the structure factor has not been clarified yet.
In this article, we combine two computational methods, molecular dynamics (MD) simulations and persistent homology (PH) in topological data analysis, to identify hidden MRO features in glassy water systems. The detected MRO structures are rings determined by the PH computing mechanism and are distinguished by their geometry, such as their size, shape and the number of components: MRO rings in LDA ice have larger size and flatter shapes, while involving less number of member atoms, than those in HDA ice. The ensemble of these MRO rings extracted from PH computations successfully reproduces the FSDP (both peak position and width) in the x-ray structure factor, demonstrating the ring structures as the real space origin of experimentally observed FSDPs. This result highlights that PH is a convincing tool to reveal a submerged order beyond the short local range in amorphous systems.

Molecular dynamics simulations
For MD simulations we used the TIP4P/2005 model [23] in the LAMMPS software [24]. We first prepared a liquid water structure as an isothermal-isobaric NPT ensemble of 2880 water molecules in a cubic box with periodic boundary conditions at T = 280 K and P = 0.1 MPa. The isobaric cooling cycle on the liquid water to obtain hyperquenched glassy water, i.e. LDA ice, ran for 200 ns at a cooling rate of 1 K ns −1 . Then it is followed by the isothermal compression of LDA ice to get HDA ice at 1.3 GPa, running for 18 ns with a compression rate of 0.1 GPa ns −1 . The time step for the integration algorithm was set to be 1 fs and the damping parameters for the thermostat and barostat are 0.5 and 5 ps, respectively. We cut off the Coulombic and Lennard-Jones interactions at distance 10 Å and 10.3092 Å , respectively [18].

Structure factor
The structure factor S(Q) measures the scattered intensity of an incident radiation for the scattering vector Q and is related to the radial distribution function (or pair correlation function) via the Fourier transform. For a perfect crystal, it exhibits infinitely many sharp peaks reflecting the periodicity of the system; for a liquid-like arrangement, a certain degree of fluctuations appears in the distribution function, yet the sharpness is weak. The numerical estimation of the structure factor of a multiatomic system of N atoms α j (1 j N) at positions x j is as follows [25]: where Q = |Q| = 2π λ is the wave number for the wavelength λ of the plane wave and is the (approximated) atomic form (or scattering) factor of each atom α j with constants a j , b j and c j reflecting the electron density of atom α j under the assumption that the electron density is spherically symmetric [26].

Topological analysis via persistent homology
For structural analysis, we used persistent homology (PH) computations from an emerging research field, topological data analysis [27,28]. PH has been actively utilized in a wide range of domains, including materials science and neuroscience [29][30][31][32][33][34], as an effective detector of topological features, i.e. holes embedded in a system. For each dimension n, we define the nth persistent homology group PH n , which extracts qualitative features of n-dimensional holes, e.g. connected components for n = 0, rings for n = 1 and cavities for n = 2, that persist over multiple scales. Since we aim to detect MRO ring structures in water that are responsible for the FSDP in the structure factor, our focus is on the case n = 1.
As one of the most efficient representations of PH computations, a persistence diagram (PD) is used to visualize the distributions of size and shape of hole structures. PD can precisely reveal hierarchical structural units in various length scales from short-to long-range order in the system. In this sense, PD has great merit compared to other conventional methods for structural analyses, such as radial distribution functions (RDF), distributions of bond or dihedral angles and ring statistics [35]. The detailed PH computing mechanism for our study using HomCloud software [36,37] follows in section 3.2 along with a schematic illustration.

Data preparation and structure factor
We performed out-of-equilibrium MD simulations on the system of 2880 water molecules to obtain LDA and HDA ices, separated by two steps (figure 1(a)) as described in section 2.1. Figure 1(b) shows the change in density at each step, which is in good agreement with the experimental densities ρ LDA ≈ 0.95 g cm −3 and ρ HDA ≈ 1.31 g cm −3 at T = 77 K and P = 1 GPa [5]. Note that we took LDA and HDA ice structures at P = 0.1 and 1300 MPa, respectively, of which densities agree the best with the experimental results.
Next, the structure factor S(Q) of each glassy water system is computed and compared with experimental results [22] to confirm the reliability of our MD-generated amorphous ice systems. (See section 2.2 for the computational details.) For an amorphous configuration, the FSDP at small Q = 0 in S(Q), if it exists, is closely related to MRO within the system [14,38]. Figure 2 compares the experimental structure factor (top) with the computed S(Q) (bottom) of our water structures generated by MD simulations. The FSDP position Q * in each case is well-matched to the experimental result at Q * LDA ≈ 1.75 and Q * HDA ≈ 2.25 within the error bound of ∼ 0.05 Å −1 . Note that the position of the peak varies slightly from the experimental result due to the difference in external conditions of temperature and pressure in our simulation setup.

Persistent homology computations
The PH computing mechanism is as follows. Suppose we have an atomic configuration of N atoms α j centered at x j ∈ R 3 (1 j N) and a set {r j ∈ R : 1 j N} of input radii. To compute PH 1 , we place a spherical ball of radius r j centered at x j , and enlarge the ball radius to r 2 j + ε with increasing ε, the scale parameter in the computation. When two balls meet, we add a line segment to connect the corresponding two atoms. As ε increases, more line segments appear and eventually create a ring by connecting multiple atoms end to end. When a ring R is born, ε is recorded as the birth scale b R of the ring. In contrast, we define the death of ring R when R is completely triangulated and every three balls has nonempty intersection. The ε at the death point is recorded as the death scale d R of the ring. The very last triangulated ring is called the death simplex, and its maximum side length represents the size of the ring. [31,39] Figure 3(a) shows a schematic illustration for the PH computing process with the resulting PD of a 2D version of a water configuration. In our study, we use the following input radii of atoms: taking the first peak positions r OH = 0.95 Å and r HH = 1.55 Å of the RDF of OH and HH pairs into consideration (figure S1 in supplementary data (stacks.iop.org/JPhysCM/31/455403/mmedia)), we solve the equations r HH = 2 r H and r OH = r O + r H to take r H = 0.775 Å and r O = 0.175 Å as the input radii of H and O atoms, respectively. As ε increases, growing balls intersect, thereby line segments connect atoms, so that multiple rings R j ( j = 1, · · · , 8) are created. Note that red triangles R j ( j = 4) are small enough to disappear shortly after their birth, while the blue ring R 4 persists longest. In the PD, the longer a ring persists, the farther away it is placed from the diagonal. The largest ring R 4 has R j for j = 5, 6, 7 as its subrings and the unfilled blue triangle at ε = ε 5 is the death simplex of R 4 . Next, the PD 1 for our MD-generated amorphous ices are aligned in figure 3(b). Note that the birth-death pairs of rings in PDs are plotted in a squared scale of Å 2 . PDs of both LDA and HDA ices show a division of points into distinct areas, characterized by lines or band-like regions, where a large number of points (each representing a ring) are concentrated. We denote the separated characteristic regions by L j and H j ( j = 1, 2, 3, 4) for LDA and HDA, respectively, as depicted in the figure. Note that rings in vertical island regions L 1 and H 1 are composed of more than four member atoms (m 4) and their triangular or quadrilateral subrings are recorded in other regions close to the diagonal. Since rings in L 4 and H 4 are small triangles whose length scales contribute little to the FSDP in the structure factor, our focus is rather on rings in other characteristic regions ( j = 1, 2, 3) that would suitably correspond to the MRO in glassy water and are likely to be associated to the sharp peaks in S(Q) in figure 2.

Characterizing MRO in LDA and HDA ices
The geometry of MRO ring structures of LDA and HDA ices substantially differs from each other in several aspects. First, stacked bar charts in figure 4(a) compare the distributions of the number of ring members (m 30) in the vertical island regions L 1 and H 1 in PDs. We denote a ring with m members by m-ring. While L 1 has a sharp distribution of the number of members, mostly concentrated on 7-10 members, H 1 has a broader distribution, as many rings of H 1 involve a larger number of atoms as their components. The median number of ring members is 9 and 12 in L 1 and H 1 , respectively. As mentioned earlier, rings in other characteristic regions for j = 1 (close to the diagonal) are simple triangles or 4-rings. We observe in figure 4(a) that HH atom pairs, as opposed to OH or OO (none) pairs, dominantly determine death scales of MRO rings both in LDA and HDA ices. The inset figures show six types of representative m-rings in L 1 ( 1 -3 ) and in H 1 ( 4 -6 ) extracted from the marked spots in the PDs in figure 3(b) with red lines connecting the atom pairs responsible for the death of each ring.
In addition to the number of ring members, MRO rings in LDA and HDA ices are differentiated by their size and shape. For this purpose, we consider the following length of ring R: which is the distance between the centers of atoms α and α , the pair of which determines the death scale d R of ring R (i.e. depicted as red lines in figure 3(a)) [39]. Note that δ(R) represents the size of R, as δ(R) is the maximum edge length within the triangulated ring, and is invariant under the choice of input radii, since it is determined by the given configuration itself. We plot red line graphs in figure 4(b) to show the mean value of δ(R) of m-rings in L 1 and H 1 regions with respect to the number of ring members m, tending to increase as m becomes larger. The average sizes δ of the rings in L 1 and H 1 contributory to the FSDPs are estimated to be notably different as 3.58 ± 0.23 Å and 2.99 ± 0.09 Å , respectively. It might be counterintuitive that rings with more members in H 1 are smaller than rings in L 1 with fewer members. This is also confirmed in PDs where the death scale (i.e. the size of the ring) exhibits a range of smaller length scales in HDA than in LDA ice. However, as water molecules become closer to each other in HDA, mainly by the pressure effect, the movement naturally leads to smaller voids in twisted ring bodies with more members, which is consistent with higher density of HDA. We define the following formula to study the degree of flatness (or planarity) of a detected MRO ring R of m members: where x j is the position of atom α j and P is the best fit plane to the ring geometry that minimizes the sum of the distances of member atoms from the plane. Note that F 0 with F = 0 indicating that all the ring members are coplanar, hence a flat ring. MRO rings in L 1 display flatter shapes with F = 0.42 ± 0.12 and F max = 0.52, while those in H 1 have F = 0.56 ± 0.19 and F max = 0.87. The mean flatness F of rings in L 1 (left) and H 1 (right) are shown with blue dotted lines in figure 4(b). Two 9-rings, each found in L 1 and H 1 regions, are shown in figure 4(c). This geometry is compatible to the picture where water molecules in the second coordination shell in the basic tetrahedral network penetrate into the first shell by the increased pressure, resulting in the deformed tetrahedral network in HDA ice.

Reproducing FSDPs by PD 1 results
To reproduce the FSDP in S(Q) from PD 1 results, we use the following equation [31]: where Π is a characteristic region L j or H j ( j = 1, 2, 3) on PDs with |Π| the number of rings in the region, (b R , d R ) is the birth-death pair of ring R in a squared scale of Å 2 and σ = 0.05 applies a smoothing effect on the curve. Figure 5 shows the reproduced FSDPs of LDA (top) and HDA (bottom) by S * (Q). We first consider each individual characteristic region Π j ( j = 1, 2, 3) on PD to plot the corresponding S * (Q) and then their combination Π 1 ∪ Π 2 ∪ Π 3 for each LDA and HDA ice. Note that rings in L 1 and L 2 regions contribute much to reproducing the peak at Q * LDA apart from the rings in L 3 region, which are of smaller scale than those in other regions. This signifies that MRO rings in L 1 and L 2 are the main contributors to the FSDP in S(Q). Also in the case of HDA, the reproduced peak from the combined characteristic region H 1 ∪ H 2 ∪ H 3 matches fairly well with the designated spot. In addition to the peak position, the width of the FSDP is also examined. We calculate the ratio Ω = w Q * HDA /w Q * LDA , where w Q * LDA and w Q * HDA are the full widths at half maximum of the FSDPs of the amorphous ices. The ratio Ω S * (Q) of the reproduced peaks is estimated to be 1.63, which excellently agrees with Ω x−ray ≈ 1.59 of the experimental data (top panel in figure 2) and Ω MD ≈ 1.60 of our MD-generated ice configurations (bottom panel in figure 2). All these values consistently reflect that the FSDP width of HDA ice is larger than that of LDA ice. Overall, PH-extracted MRO rings successfully reproduce the quantities associated with the FSDPs, including the peak positions and widths, in the structure factor, which convincingly supports that these rings construct the hidden MRO in glassy water, leading to the appearance of the sharp peaks.
In a multiatomic system as water, the larger the electron density of an atom has, the more effect it adds to the diffraction pattern. This fact has been one of the reasons for higher attention to the oxygen arrangement only in the previous studies of water structures [17,40]. As a comparison, we further carried out the FSDP reconstruction process based on the oxygen arrangement only (O-system hereafter) and compared the results to the case of the entire H 2 O system (figure S2 in supplementary data). While the trend of the FSDP of HDA appearing at larger Q than that of LDA is preserved, O-system does not accurately reproduce the quantities associated with the FSDP (e.g. peak position and width). For the O-systems in the bottom panel in figure S2, the peak positions deviate by approximately 0.4-0.6 Å −1 , and also the notably larger peak width of HDA than that of LDA is not reproduced. Thus, it is critical to keep both oxygen and hydrogen atoms into consideration for an accurate PH analysis of water in ensemble. This is a hitherto unknown result, which would have been difficult to obtain in earlier approaches, as previously known structural information was mostly based on oxygen atoms in water, overlooking hydrogen atoms as a mere decoration.

Conclusion
In conclusion, our combined approach using MD and PH has successfully embodied the real space origin of the FSDP in the structure factor of both types of amorphous ices by detecting MRO ring structures in each. We identified the MRO structures via separated regions in PDs and compared those in each type of the amorphous ice systems in terms of their shape and size as well as the number of ring constituents. The average ring size δ and the average degree of flatness F of the rings can be summarized as δ LDA = 3.58Å > 2.99Å = δ HDA and F LDA = 0.42 < 0.56 = F HDA . The PH-extracted MRO rings well explain the quantitative features (peak position and width) of the FSDPs in the structure factor, hence suitably serve as a hidden order in real space. Further applications of PH will shed light on our fundamental understanding of various types of amorphous systems with interesting MRO features that are yet hidden in a complicated geometry.