- Email: [email protected]
Two-domain structure and dynamics heterogeneity in a liquid SiO2
Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx
Contents lists available at ScienceDirect
Journal of Non-Crystalline Solids journal homepage: www.elsevier.com/locate/jnoncrysol
Two-domain structure and dynamics heterogeneity in a liquid SiO2 ⁎
P.K. Hunga, L.T. Vinha,b, , N.V. Hongc, N.T. Thu Hac, Toshiaki Iitakad a
Simulation in Materials Science Research Group, Advanced Institute of Materials Science, Ton Duc Thang University, Ho Chi Minh City, Viet Nam Faculty of Electrical and Electronics Engineering, Ton Duc Thang University, Ho Chi Minh City, Viet Nam Department of Computational Physics, Hanoi University of Science and Technology, Viet Nam d Computational Astrophysics Laboratory, RIKEN, Saitama 351-0198, Japan b c
A R T I C L E I N F O
A B S T R A C T
Keywords: Simulation Diffusion mechanism Link-cluster function Dynamics heterogeneity
Molecular dynamics simulation is employed to investigate the diffusion mechanism, dynamics and structure heterogeneity in a liquid SiO2 at ambient pressure and different temperatures. To explore these phenomena we propose the link-cluster and init-bond functions. The simulation shows that during a moderately long time the liquid has a two-domain structure that consists of separate mobile and immobile domains. These types of domains differ strongly in the rate of atomic rearrangement and fraction of defective atoms. The simulation also reveals a strong correlation between mobility of atoms and init-bond function which evidences the bondbreaking mechanism and that the reactions are non-uniformly distributed in the space. The liquid exhibits dynamics heterogeneity which is accompanied with the structure heterogeneity and cooperative movement of atoms via like-molecules SixOy.
1. Introduction Dynamics in certain liquids have been proven to be heterogeneous which concerns subsets of particles that rotate or translate farther or shorter distances than the average distance traveled by a particle with a time scale of the order of structural relaxation time. Such dynamical heterogeneity (DH) is directly related to the spatial distribution of fast and slow particles which plays a central role in several theories of glass transition. Previous studies reported that DH strongly depends on the temperature and may have consequences to a dynamical slowdown. In fact, the dramatic dynamical slowing down occurs near the glass transition point, while accompanying no noticeable change in the statistic structure [1–4]. Experimentally, DH was studied by various techniques such as the multi-dimensional nuclear resonance magnetic spectroscopy, optical photo bleaching and non-resonant spectral hole burning [5–9]. Also DH has been clarified in molecular dynamics (MD) simulations through variables such as the trajectory patterns, nonGaussian parameter, bond-lifetime patterns and multi-point dynamic susceptibilities [10–19]. The simulations in ref. [20] emerged a picture that the regions of high medium-range crystal-like order with high Debye–Waller factor form slow regions which may be an origin of DH. An opportunity to providing a potential bridge between computationally and experimentally observed DH lies in confocal microscopy experiments on dense colloidal suspensions. With the confocal microscopy, it is possible to track the motion of individual nanometer- to
⁎
micrometer-size colloidal particles at densities approaching to glass transition density, and whereby to examine the correlated motion of particles directly in the real space [21–23]. As a typical sample of network-forming liquid, silica continues drawing great interest from both practical and scientific viewpoints [24,25]. Many aspects of DH in this liquid have been successfully studied by Van Hove function and multi-point correlation functions. However the temporal–spatial distribution of fast and slow atoms accompanied with the structure heterogeneity is still poorly understood. Consequently, the picture of DH remains not fully clarified yet. Therefore the present simulation is aimed to address mentioned above problems. To do this we propose the link-cluster function characterized the clustering of atoms, and init-bond function characterized the time evolution of Si-O bond. Then we make an analysis based on those functions to give new insight into DH and structure heterogeneity. Furthermore we focus on the bond-breaking mechanism occurred in the liquid silica. 2. Computational procedure We have prepared models consisting of 1000 Si and 2000 O by means of MD simulation. The van Beest–Kramer–van Santen (BKS) potential is employed [26]. Initial configuration is generated by randomly placing all atoms in a simple box. The obtained sample is heated to 5000 K and then cooled down to specified temperatures at ambient
Corresponding author at: Faculty of Electrical and Electronics Engineering, Ton Duc Thang University, Ho Chi Minh City, Viet Nam. E-mail addresses: [email protected] (P.K. Hung), [email protected] (L.T. Vinh).
https://doi.org/10.1016/j.jnoncrysol.2018.01.023 Received 17 November 2017; Received in revised form 10 January 2018; Accepted 15 January 2018 0022-3093/ © 2018 Elsevier B.V. All rights reserved.
Please cite this article as: Hung, P.K., Journal of Non-Crystalline Solids (2018), https://doi.org/10.1016/j.jnoncrysol.2018.01.023
Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx
P.K. Hung et al.
pressure. After that the sample has been relaxed in NVE ensemble (the constant volume and energy) until reach the equilibrium. Two samples have been prepared at temperature of 3000 and 3500 K. The structural characteristics of these samples are in good agreement with experiment. More details about preparing of silica model can be found elsewhere [27,28]. We also prepared a large model of 4998 atoms (1666 Si and 3332 O) at 3000 K and ambient pressure. The calculation shows no noticeable difference in dynamics and structure between 3000-atom and 4998-atom models. Consider a set of N atoms during a time tobs. The temporal-spatial distribution of these atoms was analyzed by link-cluster functions. We use following notations and definitions. The linkages between atoms are determined by a distance rlk. Two atoms form a linkage if the distance between them is smaller than rlk. The link-cluster is a subset of atoms connected to each other through linkages. It means that for any two atoms k1 and ks of the cluster one can find other atoms k2, …, ks−1 of this cluster so that the distance between atoms ki and ki+1 is smaller than rlk where i = 1, 2, … s − 1. The link-cluster function Flink(r, t) is defined as a number of link-clusters formed by atoms of the given set at a time t and rlk = r, where t ≤ tobs. For the set of N atoms Flink(r, t) varies from 1 to N. The calculation of Flink(r, t) is performed as follows. Firstly, we find linkages for all atoms. Then each atom has been preliminarily assigned to a label k, where k = 1, 2, … N. Next, the atoms k1 and k2 are reassigned to a new label, if they form a linkage and k1 ≠ k2. Both atoms are reassigned to k1 if k1 < k2, otherwise to k2. This procedure is performed until all pairs of atoms forming a linkage have the same label. Finally, we obtain a list of labels [k1, k2, …km] and link-clusters [Cl1, Cl2, …, Clm], where Flink(r, t) = m; m is the number of link-cluster. The link-cluster Cli consists of atoms having the label ki. We consider sets of mobile, immobile and random atom (SMA, SIMA and SRA). Each set contains 600 atoms, i.e. 20% of total atoms. SMA comprises atoms that have qit2 larger than remaining atoms, where qit2 = [qi(0) − qi(t)]2 is the square displacement of ith atom; qi(t) is the coordinate of ith atom at the time t. In contrast, SIMA consists of atoms of which qit2 is smaller than remaining atoms. Here all atoms of the sample are divided in two groups: SMA (or SIMA) and remaining atoms. The atoms of SRA are randomly chosen from the sample. For the simplicity we call mobile and immobile atoms as selective atoms. The list of selective and random atoms is determined by a MD observation time tobs. We examine trajectories of atoms of this list within an interval [0, tobs]. For a moderately long tobs the spatial distributions are different for selective and random atoms. This difference demonstrates DH in the liquid. As tobs increases, the list of atoms from considered sets changes, and as we show below, the temporal-spatial distribution of selective atoms approaches to the one of random atoms. As shown in refs. [24, 25], the structure of silica liquid consists of a network of SiO4 connected by bridging oxygen. There is also a small amount of defective units such as SiO5. We call this network as Si-O network. For the samples at 3000 and 3500 K, the fraction of SiO4 is found to be 0.92 and 0.83, respectively. The fraction of OSi2 for those samples is equal to 0.94 and 0.90. The diffusion follows the bondbreaking mechanism [27, 28]. The breaking and reforming of SieO bond are realized via reactions:
SiO4 ↔ SiO5 and OSi2 ↔ OSi3
Fig. 1. The link-cluster function for two moments: t = 71.7 (a) and 143.4 ps (b). The time tobs is equal to 143.4 ps. It can be seen that Flink(r,t) drops drastically as r varies from 1.3 to 1.9 Å.
(1)
The cutoff distance rcutoff is set to 2.45 Å which is chosen from first minimum of the pair radial distribution function. If the distance between O and Si is less than rcutoff, then they form a SieO bond. The Si-O bonds in the starting configuration (t = 0) are called init-bonds. The number of init-bond for a configuration at t > 0 is expressed as
Finb (t ) =
Fig. 2. The number of SieO bonds as a function of time. Here SMA-SIMA curve corresponds to the bond formed by one atom belonging to SMA and another one of SIMA.
∑ δ (rijo , rijt) i
3. Results and discussion (2) Fig. 1 displays Flink(r, t). The time tobs equals 143.4 ps. At 3000 K and as r varies from 1.3 to 1.9 Å, Flink(r, t) for SIMA, SMA and SRA drops drastically to 250, 390 and 450, respectively. With further increasing r
where rij0, rijt is the distance between ith and jth atoms of configurations at t = 0 and t > 0, respectively; δ (rij0,rijt) = 1 if rij0 < rcutoff and rijt < rcutoff, otherwise δ (rij0,rijt) = 0. 2
Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx
P.K. Hung et al.
a shoulder appears, and then Flink(r, t) decreases gradually. At 3500 K Flink(r, t) is similar for SMA and SRA. Moreover Flink(r, t) for SRA differs from the one for SIMA, but the difference is less pronounced compared to the case of 3000 K. We conclude that at low temperature the selective atoms are non-uniformly distributed in the space, while the spatial distribution of random atoms is more homogeneous. As shown from Fig. 1, Flink(r, t) is almost unchanged for the cases t = 71.7 and 143.4 ps, although the atoms of considered sets move in space and time. Therefore, the temporal-spatial distribution of selective and random atoms slightly changes during the time tobs. We calculate numbers of SieO bonds formed by atoms from considered sets. The result is presented in Fig. 2. At 3000 K the curves are completely separate and the number of bonds decreases in the order: SIMA → SMA → SRA. In the case of 3500 K the number of bonds is close for SMA and SIMA, but the one for SIMA remains largest. Thus, the selective atoms tend to concentrate in some domains to form SieO bonds. Consequently, the SieO network contains large SixOy subnets of selective atoms which are schematically illustrated in Fig. 3. Because of these subnets belong to link-clusters with rlk = 1.9 Å, the observed drop of Flink(r, t) is due to SixOy subnets. The shoulder seen in Fig. 1 is related to that the number of pairs of atoms separated at a distance from 1.9 to 2.6 Å is relatively small. This can be seen from pair radial distribution function shown in Fig. 4. Accordingly, the number of SieO pairs separated from 1.9 to 2.6 Å is significantly smaller than the one from 1.3 to 1.9 Å. In Fig. 2 we also plot the number of SieO bonds formed by two atoms: one belonging to SMA and another to SIMA (SMA-SIMA type). We note that at 3000 K the number of bonds of this type is very small. This indicates that most of mobile atoms are located outside domains in which the immobile atoms are concentrated. The behavior of selective atoms demonstrates that the SieO network contains separate domains in which the atoms rearrange much weaker than in a remaining part of the network. We call those domains immobile domains. The remaining part is called a mobile domain. Unlike the immobile domain, the atoms make a strong rearrangement in the mobile domain. It turns out that the liquid has two-domain structure which is schematically illustrated in Fig. 5. For the convenience of discussion we use following notations. The density of atoms from considered sets is denoted by CCS. Here CCS = 600 / 3000 = 0.2. The average density of immobile atoms in immobile domains equals CIMD = < njimd / ntjimd > where njimd, ntjimd is the number of immobile atoms and total atoms in jth immobile domain,
Fig. 3. Schematics of SixOy subnets of selective and random atoms; the black line represents the SieO bond; the red and blue circle represents Si and O, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 4. The pair radial distribution function for SieO pair.
Fig. 5. The two-domain structure of the liquid for the case of low (a) and high temperature (b). The liquid contains a number of immobile domains in which the atoms rearrange weakly. In contrast, the atomic rearrangement is strong in a mobile domain.
3
Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx
P.K. Hung et al.
respectively. The brackets mean averaging over all immobile domains. We denote the density of mobile atoms in the mobile domain by CMD. This density is equal to nmd/ntmd where nmd, ntmd is the number of mobile atoms and total atoms in the mobile domain, respectively. Because of immobile atoms are concentrated in immobile domains and most of mobile atoms reside in the mobile domain, we have: CIMD > CCS and CMD > CCS. This leads to that large SixOy subnets of immobile and mobile atom are formed in the separate domains as observed from the simulation. Moreover, due to CIMD > CMD the number of bonds for immobile atoms is larger than mobile atoms. With increasing temperature, the size of immobile domains decreases and CMD approaches to CCS. To get more details about spatial distributions of selective atoms we calculate the size distribution of link-clusters with rlk = 1.9 Å. The size of link-cluster is defined as a number of atoms belonging to this cluster. The result is presented in Tables 1 and 2. At 3000 K the size of linkcluster for SRA varies from 1 to 6, whereas SMA and SIMA have numerous clusters with size from 7 to 36. The average size of link-cluster for SIMA is significantly larger than that for SMA. At 3500 K those distributions are similar for SRA and SMA which means that mobile atoms are uniformly distributed in the space at high temperature. Moreover they show a clear difference between SMA and SIMA. However, the difference is less pronounced compared to the case of 3000 K. The result seen in Tables 1 and 2 can be interpreted by that immobile atoms are concentrated in immobile domains which leads to forming large link-clusters of immobile atom. In contrast, the majority of mobile atoms reside in the mobile domain so that large link-clusters of mobile atom are present in the mobile domain. The spatial distribution of selective and random atoms can be described schematically in Fig. 6. It can be seen that there are numerous large link-clusters of selective atoms, while random atoms are distributed in the space by small clusters. The system contains only small clusters of random atoms because of CIMD > CCS and CMD > CCS. The average size of link-cluster of immobile atom is significantly larger than the one of mobile atom due to CIMD > CMD. The size distribution of link-cluster is almost unchanged with time because of CIMD and CMD are slightly fluctuated during the time tobs. At 3500 K the size distribution of link-cluster is similar for SRA and SMA due to that the average size of immobile
Table 1 The number of link-clusters in configurations at t = 71.7 and 143.4 ps, and at 3000 K. Here SCl, NCl is the size of link-cluster and number of clusters having the corresponding SCl, respectively. SCl
1 2 3 4 5 6 7 8 9 11 16 17 19 20 36 –
NCl, t = 71.7 ps
SCl
SMA
SIMA
SRA
310 14 22 10 8 4 3 3 1 2 1 0 0 0 0 –
135 61 18 8 6 5 2 8 3 0 0 1 1 1 1 –
338 75 22 6 2 2 0 0 0 0 0 0 0 0 0 –
1 2 3 4 5 6 7 8 9 10 11 13 14 16 23 30
NCl, t = 143.4 ps SMA
SIMA
SRA
312 27 23 9 2 4 2 2 1 2 0 1 0 0 1 0
149 49 17 10 9 5 2 8 2 2 1 0 1 1 0 1
345 68 26 5 3 1 0 0 0 0 0 0 0 0 0 0
Table 2 The number of link-clusters in the configuration at t = 71.7 and 143.4 ps, and at 3500 K. SCl
1 2 3 4 5 6 7 8
NCl, t = 71.7 ps
SCl
SMA
SIMA
SRA
396 30 24 15 1 0 1 0
268 59 28 16 5 2 3 1
335 76 22 6 2 1 1 0
1 2 3 4 5 6 7 8
NCl, t = 143.4 ps SMA
SIMA
SRA
395 37 20 5 4 4 1 0
307 65 14 9 6 4 1 3
361 61 29 5 2 0 0 0
Fig. 6. The spatial distribution of random and selective atoms. Here we represent the atom by a circle with radius of rlk / 2. If two atoms form a linkage, then the circles overlap with each other. A linkcluster therefore corresponds to a subset of overlapped circles. The red, yellow and blue circle represents the random, immobile and mobile atom, respectively. The other-type atoms are not shown here. It can be seen that selective atoms form numerous large linkclusters, whereas random atoms form only small ones. The average size of link-cluster of immobile atoms is larger than that of mobile atoms. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
4
Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx
P.K. Hung et al.
Table 3 The time evolution of link-cluster with size larger than four for SMA at 3000 K. t=0
t = 14.34 ps
t = 71.7 ps
t = 143.4 ps
SCl
NCl
SCl
NCl
SCl
NCl
SCl
NCl
5 6 7 8 9 12 15 – – – – –
6 1 2 3 2 1 1 – – – – –
1 3 4 5 6 7 8 9 12 16 – –
3 1 3 4 2 3 2 2 1 1 – –
1 2 3 4 5 6 7 10 11 12 – –
15 5 7 7 3 2 3 1 2 2 – –
1 2 3 4 5 6 7 8 9 10 13 23
24 9 8 4 2 3 1 1 1 2 1 1
Table 4 The time evolution of link-cluster with size larger than four for SIMA at 3000 K. t=0
t = 14.34 ps
t = 71.7 ps
Fig. 7. The number of defective atoms as a function of time.
t = 143.4 ps
SCl
NCl
SCl
NCl
SCl
NCl
SCl
NCl
5 6 7 8 9 11 17 18 – – – – – – –
7 7 2 7 3 1 2 2 – – – – – – –
1 2 3 4 5 6 7 8 9 10 11 12 17 18 20
3 1 1 3 6 5 1 7 4 1 1 1 1 1 1
1 2 3 5 6 7 8 9 10 14 17 19 34 – –
3 1 1 6 3 2 9 2 1 1 1 2 1 – –
1 2 3 4 5 6 7 8 9 10 11 14 16 30 –
6 4 1 2 6 5 2 8 2 2 1 1 1 1 –
of the found list and may consist of other selective atoms. Finally, we determine reforming clusters in different configurations at t > 0. The calculation result for four configurations is shown in Tables 3 and 4. We note that there are a number of small reforming clusters with size less than four. For SMA the number of small cluster largely increases with time (see Table 5). This indicates that initial clusters strongly split into small clusters. Therefore the atoms make a strong rearrangement in the mobile domain. On other hands, the number of large clusters with size bigger than four slightly decreases from 16 to 12. The number of total atoms belonging to large clusters also weakly changes. In particular, NTA varies from 104 to 115 that slightly smaller than the one of initial clusters equal to 119. From this follows that initial clusters are split into small clusters which then merge with other atoms to form large clusters with size up to 23 atoms. This evidences that CMD is much larger than CCS. It is expected that the mobile atoms are redistributed from mobile to immobile domains due to CIMD > CMD. Consequently, CMD and the average size of reforming clusters decreases as t increases. However, the average size of reforming clusters NTA/NLCl slightly varies during the time tobs. This evidences that the exchange of mobile atoms between mobile and immobile domains proceeds so that CMD is slightly fluctuated. In the case of SIMA the number of small cluster increases from 8 to 13 that significantly smaller compared to the case of SMA. There are numerous large clusters with size up to 34 atoms. NTA/NLCl is also slightly changes with time. This result evidences CIMD slightly varies during the time tobs. Also the rate of atomic rearrangement in the immobile domain is significantly smaller than the one in the mobile domain. We calculate the number of defective atoms from considered sets. The defective atom is defined as over-coordinated Si or O. Most of defective atoms have the coordination number of five for Si and three for O. The Si with the coordination number of four and O having two coordinated atoms is called normal atom. In Fig. 7 we plot the number of defective atom as a function of time. At 3000 K the number of defective atoms of SMA is significantly larger than the one of SIMA. This is originated from that the diffusion is mainly realized by reactions (1) via defective units such as SiO5 and OSi3. We conclude that the immobile domains have large fraction of normal atom. In contrast, the mobile domain possesses large fraction of defective atoms. Therefore the immobile and mobile domains possess different local microstructures
domains significantly decreases and CMD is close to CCS. This result indicates that the liquid exhibits DH the degree of which reduces with increasing temperature. To get more information about mobile and immobile domains we examine the evolution of large link-clusters with rlk = 1.9 Å. The calculation is performed as follows. Firstly, we determine all link-clusters with size larger than four in the starting configuration (t = 0). The clusters found are called initial clusters. Next, a list of atoms belonging to initial clusters is determined. As the time proceeds, the atoms of the found list are redistributed through new link-clusters called reforming clusters. Unlike initial cluster, the reforming cluster has at least an atom Table 5 Characteristics of initial and reforming clusters. Here NSCl, NLCl is the number of cluster with size SCl < 5 and SCl ≥ 5, respectively; NTA is the number of total atoms belonging to the clusters with size SCl ≥ 5. t, ps
0 14.34 71.7 143.4
SMA
SIMA
NSCl
NLCl
NTA
NSCl
NLCl
NTA
0 7 34 45
16 15 13 12
119 115 104 108
0 8 11 13
31 29 28 29
255 247 265 247
5
Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx
P.K. Hung et al.
Fig. 8. The link-cluster function of SMA and SIMA for tobs = 71.7 and 143.4 ps at 3000 (a) and 3500 K (b).
two-domain structure is observed within a time interval [0, tobs]. For the next interval [tobs, 2tobs] we observe other immobile domains which possess similar properties. In particular, the size distribution of immobile domain is almost unchanged, but their shape and position change significantly. Overall, the scenario for dynamics can be described as follows. During a moderate time tobs the immobile atoms are concentrated in immobile domains, which lead to forming large link-clusters with rlk = 1.9 Å. Meanwhile most of mobile atoms move inside a mobile domain so that large link-clusters of mobile atom are formed there. The rate of atomic rearrangement in the immobile domain is significantly smaller than the one in the mobile domain. The exchange of selective
which mean that DH is accompanied with the structure heterogeneity in the liquid. In the case of 3500 K the number of defective atoms is close for SIMA and SMA. This result can be interpreted by that at this temperature the rate of reaction (1) is high and the spatial distribution of selective atoms is more homogeneous. In Fig. 8 we present Flink(r, t) for different tobs. One can see that Flink(r, t) for selective atoms approaches to the one for random atoms as the time tobs increases. This means that the spatial correlation between selective atoms is weaker as tobs is longer. As a result, the size of immobile domains decreases with increasing tobs, i.e. upon long time tobs we do not observe the two-domain structure. It is worth to note that the
6
Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx
P.K. Hung et al.
the spatial distribution of selective atoms is non-uniform, the reactions are also non-uniformly distributed in the space. Due to the SieO bond the atoms do not make random displacements, but they tend to move by a group as a like-molecule SixOy, where x > 0, y > 0. For instance, Si2O7 kept eight init-bonds. The SixOy molecule with lifetime t is defined as a subset of atoms connected to each other via init-bonds over whole time t. Table 6 presents the number of such molecules for different lifetimes. As expected, the size of large SixOy decreases as the time t increases. We observed that there are numerous large SixOy molecules existed for long time even at high temperature. This clearly demonstrates the cooperative movement of atoms via SixOy. Such movement results in that the selective atoms tend to locate nearby to form Si-O bonds. Moreover the selective atoms form numerous large link-clusters with rlk = 1.9 Å. We conclude that DH observed is accompanied with the cooperative movement via SixOy. 4. Conclusions MD simulation is carried out to study the diffusion mechanism, dynamics and structure heterogeneity in the liquid silica. We propose link-cluster and init-bond functions to explore these phenomena. The calculation shows that at low temperature selective atoms are nonuniformly distributed in the space. In particular, the selective atoms form numerous large link-clusters with rlk = 1.9 Å. At 3000 K the immobile and mobile atoms form respectively large clusters with size up to 34 and 23 atoms. Whereas the random atoms are arranged by small clusters with size < 6 atoms. We reveal that during a moderately long time tobs the liquid has a two-domain structure that consists of separate immobile and mobile domains. Immobile domains strongly differ from the mobile domain in the rate of atomic rearrangement and that the fraction of defective atoms in immobile domains is significantly larger than the one in the mobile domain. The calculation also reveals that the diffusion is realized via bondbreaking mechanism. The reactions are non-uniformly distributed in the space and the liquid exhibits DH the degree of which reduces with increasing temperature. Furthermore the simulation shows a strong correlation between mobility of atoms and init-bond function. We find that DH is accompanied with structure heterogeneity and cooperative movement of atoms via like-molecules SixOy.
Fig. 9. The Finb(r) as a function of time for SMA, SIMA and SRA.
Table 6 The number of SixOy like-molecules. Smolecule, Nmolecule is the range of molecule size and number of molecules having the corresponding sizes, respectively; T is the temperature. t = 71.7 ps; T = 3000 K
t = 143.4 ps; T = 3000 K
t = 71.7 ps; T = 3500 K
t = 143.4 ps; T = 3500 K
Smolecule
Nmolecule
Smolecule
Nmolecule
Smolecule
Nmolecule
Smolecule
Nmolecule
2–9 10–12 2469
91 3 1
2–9 10–37 1549
248 11 1
2–9 10–14 20
483 3 1
2 3 4–5
149 25 4
Acknowledgement
atoms between mobile and immobile domains proceeds so that CIMD and CMD are slightly fluctuated within the time tobs. We calculate Finb(t) by averaging over all O and Si which belong to SIMA, SMA and SRA. The result is plotted in Fig. 9. Clearly that the mobility of atom is strongly correlated with Finb(t). Namely, Finb(t) decreases in the order: SIMA → SMA → SRA. Also at 3500 K the decrease of Finb(t) with time proceeds much faster than at 3000 K. These results support the bond-breaking mechanism by the fact that the atoms possessing high mobility undergo more reactions, i.e. more init-bonds are replaced. The fast decrease of Finb(t) with time at 3500 K is caused by that the SieO bond is broken more easy at high temperature. Because
The authors are grateful for support by the NAFOSTED Vietnam (grant 103.01-2015.12). This research was supported in part by MEXT as “Exploratory Challenge on Post-K computer” (Frontiers of Basic Science: Challenging the Limits) and by RIKEN as the iTHES Project. References [1] H.C. Andersen, Proc. Natl. Acad. Sci. U. S. A. 102 (2005) 6686. [2] F. Sciortino, P. Tartaglia, Phys. 54 (2005) 471.
7
Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx
P.K. Hung et al.
[3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
[15] Sharon C. Glotzer, J. Non-Cryst. Solids 274 (2000) 342. [16] T. Kawasaki, T. Araki, H. Tanaka, Phys. Rev. Lett. 99 (2007) 215701. [17] Wen-Sheng Xu, Zhao-Yan Sun, Li-Jia An, J. Phys. Condens. Matter 24 (2012) 325101. [18] G.N. Greaves, S. Sen, Adv. Phys. 56 (2007) 1. [19] J.S. Langer, Rep. Prog. Phys. 77 (2014) 042501. [20] T. Kawasaki, H. Tanaka, J. Phys. Condens. Matter 23 (2011) 194121. [21] E.R. Weeks, J.C. Crocker, D.A. Weitz, J. Phys. Condens. Matter 19 (2007) 205131. [22] E.R. Weeks, D.A. Weitz, Phys. Rev. Lett. 89 (2002) 095704. [23] C. Donati, S.C. Glotzer, P.H. Poole, Phys. Rev. Lett. 82 (1999) 5064. [24] D. Coslovich, G. Pastore, J. Phys. Condens. Matter 21 (2009) 285107. [25] Y. Wang, et al., Nat. Commun. 5 (2014) 3241. [26] B. van Beest, G. Kramer, R. van Santen, Phys. Rev. Lett. 64 (1990) 1955. [27] P.K. Hung, L.T. Vinh, J. Non-Cryst. Solids 352 (2006) 5531. [28] P.K. Hung, N.T.T. Ha, N.V. Hong, Eur. Phys. J. E 36 (2013) 60.
R. Richert, J. Phys. Condens. Matter 14 (2002) 703. T. Kawasaki, H. Tanaka, J. Phys. Condens. Matter 22 (2010) 232102. I. Chang, H. Sillescu, J. Phys. Chem. B 101 (1997) 8794. H. Sillescu, J. Non-Cryst. Solids 243 (1999) 81. M.T. Cicerone, M.D. Ediger, J. Chem. Phys. 104 (1996) 7210. B. Schiener, R.V. Chamberlin, G. Diezemann, R. Bohmer, J. Chem. Phys. 107 (1997) 7746. R. Richert, J. Phys. Chem. B 101 (1997) 6323. M. Shajahan, G. Razul, Gurpreet S. Matharoo, Peter H. Poole, J. Phys. Condens. Matter 23 (2011) 235103. Richard K. Darst, David R. Reichman, Giulio Biroli, J. Chem. Phys. 132 (2010) 044510. Z. Zheng, et al., Nat. Commun. 5 (2014) 3829. G.A. Appignanesi, J.A. Rodriguez Fris, J. Phys. Condens. Matter 21 (2009) 203103. Tanaka, et al., Nat. Mater. 9 (2010) 324.
8
Contents lists available at ScienceDirect
Journal of Non-Crystalline Solids journal homepage: www.elsevier.com/locate/jnoncrysol
Two-domain structure and dynamics heterogeneity in a liquid SiO2 ⁎
P.K. Hunga, L.T. Vinha,b, , N.V. Hongc, N.T. Thu Hac, Toshiaki Iitakad a
Simulation in Materials Science Research Group, Advanced Institute of Materials Science, Ton Duc Thang University, Ho Chi Minh City, Viet Nam Faculty of Electrical and Electronics Engineering, Ton Duc Thang University, Ho Chi Minh City, Viet Nam Department of Computational Physics, Hanoi University of Science and Technology, Viet Nam d Computational Astrophysics Laboratory, RIKEN, Saitama 351-0198, Japan b c
A R T I C L E I N F O
A B S T R A C T
Keywords: Simulation Diffusion mechanism Link-cluster function Dynamics heterogeneity
Molecular dynamics simulation is employed to investigate the diffusion mechanism, dynamics and structure heterogeneity in a liquid SiO2 at ambient pressure and different temperatures. To explore these phenomena we propose the link-cluster and init-bond functions. The simulation shows that during a moderately long time the liquid has a two-domain structure that consists of separate mobile and immobile domains. These types of domains differ strongly in the rate of atomic rearrangement and fraction of defective atoms. The simulation also reveals a strong correlation between mobility of atoms and init-bond function which evidences the bondbreaking mechanism and that the reactions are non-uniformly distributed in the space. The liquid exhibits dynamics heterogeneity which is accompanied with the structure heterogeneity and cooperative movement of atoms via like-molecules SixOy.
1. Introduction Dynamics in certain liquids have been proven to be heterogeneous which concerns subsets of particles that rotate or translate farther or shorter distances than the average distance traveled by a particle with a time scale of the order of structural relaxation time. Such dynamical heterogeneity (DH) is directly related to the spatial distribution of fast and slow particles which plays a central role in several theories of glass transition. Previous studies reported that DH strongly depends on the temperature and may have consequences to a dynamical slowdown. In fact, the dramatic dynamical slowing down occurs near the glass transition point, while accompanying no noticeable change in the statistic structure [1–4]. Experimentally, DH was studied by various techniques such as the multi-dimensional nuclear resonance magnetic spectroscopy, optical photo bleaching and non-resonant spectral hole burning [5–9]. Also DH has been clarified in molecular dynamics (MD) simulations through variables such as the trajectory patterns, nonGaussian parameter, bond-lifetime patterns and multi-point dynamic susceptibilities [10–19]. The simulations in ref. [20] emerged a picture that the regions of high medium-range crystal-like order with high Debye–Waller factor form slow regions which may be an origin of DH. An opportunity to providing a potential bridge between computationally and experimentally observed DH lies in confocal microscopy experiments on dense colloidal suspensions. With the confocal microscopy, it is possible to track the motion of individual nanometer- to
⁎
micrometer-size colloidal particles at densities approaching to glass transition density, and whereby to examine the correlated motion of particles directly in the real space [21–23]. As a typical sample of network-forming liquid, silica continues drawing great interest from both practical and scientific viewpoints [24,25]. Many aspects of DH in this liquid have been successfully studied by Van Hove function and multi-point correlation functions. However the temporal–spatial distribution of fast and slow atoms accompanied with the structure heterogeneity is still poorly understood. Consequently, the picture of DH remains not fully clarified yet. Therefore the present simulation is aimed to address mentioned above problems. To do this we propose the link-cluster function characterized the clustering of atoms, and init-bond function characterized the time evolution of Si-O bond. Then we make an analysis based on those functions to give new insight into DH and structure heterogeneity. Furthermore we focus on the bond-breaking mechanism occurred in the liquid silica. 2. Computational procedure We have prepared models consisting of 1000 Si and 2000 O by means of MD simulation. The van Beest–Kramer–van Santen (BKS) potential is employed [26]. Initial configuration is generated by randomly placing all atoms in a simple box. The obtained sample is heated to 5000 K and then cooled down to specified temperatures at ambient
Corresponding author at: Faculty of Electrical and Electronics Engineering, Ton Duc Thang University, Ho Chi Minh City, Viet Nam. E-mail addresses: [email protected] (P.K. Hung), [email protected] (L.T. Vinh).
https://doi.org/10.1016/j.jnoncrysol.2018.01.023 Received 17 November 2017; Received in revised form 10 January 2018; Accepted 15 January 2018 0022-3093/ © 2018 Elsevier B.V. All rights reserved.
Please cite this article as: Hung, P.K., Journal of Non-Crystalline Solids (2018), https://doi.org/10.1016/j.jnoncrysol.2018.01.023
Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx
P.K. Hung et al.
pressure. After that the sample has been relaxed in NVE ensemble (the constant volume and energy) until reach the equilibrium. Two samples have been prepared at temperature of 3000 and 3500 K. The structural characteristics of these samples are in good agreement with experiment. More details about preparing of silica model can be found elsewhere [27,28]. We also prepared a large model of 4998 atoms (1666 Si and 3332 O) at 3000 K and ambient pressure. The calculation shows no noticeable difference in dynamics and structure between 3000-atom and 4998-atom models. Consider a set of N atoms during a time tobs. The temporal-spatial distribution of these atoms was analyzed by link-cluster functions. We use following notations and definitions. The linkages between atoms are determined by a distance rlk. Two atoms form a linkage if the distance between them is smaller than rlk. The link-cluster is a subset of atoms connected to each other through linkages. It means that for any two atoms k1 and ks of the cluster one can find other atoms k2, …, ks−1 of this cluster so that the distance between atoms ki and ki+1 is smaller than rlk where i = 1, 2, … s − 1. The link-cluster function Flink(r, t) is defined as a number of link-clusters formed by atoms of the given set at a time t and rlk = r, where t ≤ tobs. For the set of N atoms Flink(r, t) varies from 1 to N. The calculation of Flink(r, t) is performed as follows. Firstly, we find linkages for all atoms. Then each atom has been preliminarily assigned to a label k, where k = 1, 2, … N. Next, the atoms k1 and k2 are reassigned to a new label, if they form a linkage and k1 ≠ k2. Both atoms are reassigned to k1 if k1 < k2, otherwise to k2. This procedure is performed until all pairs of atoms forming a linkage have the same label. Finally, we obtain a list of labels [k1, k2, …km] and link-clusters [Cl1, Cl2, …, Clm], where Flink(r, t) = m; m is the number of link-cluster. The link-cluster Cli consists of atoms having the label ki. We consider sets of mobile, immobile and random atom (SMA, SIMA and SRA). Each set contains 600 atoms, i.e. 20% of total atoms. SMA comprises atoms that have qit2 larger than remaining atoms, where qit2 = [qi(0) − qi(t)]2 is the square displacement of ith atom; qi(t) is the coordinate of ith atom at the time t. In contrast, SIMA consists of atoms of which qit2 is smaller than remaining atoms. Here all atoms of the sample are divided in two groups: SMA (or SIMA) and remaining atoms. The atoms of SRA are randomly chosen from the sample. For the simplicity we call mobile and immobile atoms as selective atoms. The list of selective and random atoms is determined by a MD observation time tobs. We examine trajectories of atoms of this list within an interval [0, tobs]. For a moderately long tobs the spatial distributions are different for selective and random atoms. This difference demonstrates DH in the liquid. As tobs increases, the list of atoms from considered sets changes, and as we show below, the temporal-spatial distribution of selective atoms approaches to the one of random atoms. As shown in refs. [24, 25], the structure of silica liquid consists of a network of SiO4 connected by bridging oxygen. There is also a small amount of defective units such as SiO5. We call this network as Si-O network. For the samples at 3000 and 3500 K, the fraction of SiO4 is found to be 0.92 and 0.83, respectively. The fraction of OSi2 for those samples is equal to 0.94 and 0.90. The diffusion follows the bondbreaking mechanism [27, 28]. The breaking and reforming of SieO bond are realized via reactions:
SiO4 ↔ SiO5 and OSi2 ↔ OSi3
Fig. 1. The link-cluster function for two moments: t = 71.7 (a) and 143.4 ps (b). The time tobs is equal to 143.4 ps. It can be seen that Flink(r,t) drops drastically as r varies from 1.3 to 1.9 Å.
(1)
The cutoff distance rcutoff is set to 2.45 Å which is chosen from first minimum of the pair radial distribution function. If the distance between O and Si is less than rcutoff, then they form a SieO bond. The Si-O bonds in the starting configuration (t = 0) are called init-bonds. The number of init-bond for a configuration at t > 0 is expressed as
Finb (t ) =
Fig. 2. The number of SieO bonds as a function of time. Here SMA-SIMA curve corresponds to the bond formed by one atom belonging to SMA and another one of SIMA.
∑ δ (rijo , rijt) i
3. Results and discussion (2) Fig. 1 displays Flink(r, t). The time tobs equals 143.4 ps. At 3000 K and as r varies from 1.3 to 1.9 Å, Flink(r, t) for SIMA, SMA and SRA drops drastically to 250, 390 and 450, respectively. With further increasing r
where rij0, rijt is the distance between ith and jth atoms of configurations at t = 0 and t > 0, respectively; δ (rij0,rijt) = 1 if rij0 < rcutoff and rijt < rcutoff, otherwise δ (rij0,rijt) = 0. 2
Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx
P.K. Hung et al.
a shoulder appears, and then Flink(r, t) decreases gradually. At 3500 K Flink(r, t) is similar for SMA and SRA. Moreover Flink(r, t) for SRA differs from the one for SIMA, but the difference is less pronounced compared to the case of 3000 K. We conclude that at low temperature the selective atoms are non-uniformly distributed in the space, while the spatial distribution of random atoms is more homogeneous. As shown from Fig. 1, Flink(r, t) is almost unchanged for the cases t = 71.7 and 143.4 ps, although the atoms of considered sets move in space and time. Therefore, the temporal-spatial distribution of selective and random atoms slightly changes during the time tobs. We calculate numbers of SieO bonds formed by atoms from considered sets. The result is presented in Fig. 2. At 3000 K the curves are completely separate and the number of bonds decreases in the order: SIMA → SMA → SRA. In the case of 3500 K the number of bonds is close for SMA and SIMA, but the one for SIMA remains largest. Thus, the selective atoms tend to concentrate in some domains to form SieO bonds. Consequently, the SieO network contains large SixOy subnets of selective atoms which are schematically illustrated in Fig. 3. Because of these subnets belong to link-clusters with rlk = 1.9 Å, the observed drop of Flink(r, t) is due to SixOy subnets. The shoulder seen in Fig. 1 is related to that the number of pairs of atoms separated at a distance from 1.9 to 2.6 Å is relatively small. This can be seen from pair radial distribution function shown in Fig. 4. Accordingly, the number of SieO pairs separated from 1.9 to 2.6 Å is significantly smaller than the one from 1.3 to 1.9 Å. In Fig. 2 we also plot the number of SieO bonds formed by two atoms: one belonging to SMA and another to SIMA (SMA-SIMA type). We note that at 3000 K the number of bonds of this type is very small. This indicates that most of mobile atoms are located outside domains in which the immobile atoms are concentrated. The behavior of selective atoms demonstrates that the SieO network contains separate domains in which the atoms rearrange much weaker than in a remaining part of the network. We call those domains immobile domains. The remaining part is called a mobile domain. Unlike the immobile domain, the atoms make a strong rearrangement in the mobile domain. It turns out that the liquid has two-domain structure which is schematically illustrated in Fig. 5. For the convenience of discussion we use following notations. The density of atoms from considered sets is denoted by CCS. Here CCS = 600 / 3000 = 0.2. The average density of immobile atoms in immobile domains equals CIMD = < njimd / ntjimd > where njimd, ntjimd is the number of immobile atoms and total atoms in jth immobile domain,
Fig. 3. Schematics of SixOy subnets of selective and random atoms; the black line represents the SieO bond; the red and blue circle represents Si and O, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 4. The pair radial distribution function for SieO pair.
Fig. 5. The two-domain structure of the liquid for the case of low (a) and high temperature (b). The liquid contains a number of immobile domains in which the atoms rearrange weakly. In contrast, the atomic rearrangement is strong in a mobile domain.
3
Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx
P.K. Hung et al.
respectively. The brackets mean averaging over all immobile domains. We denote the density of mobile atoms in the mobile domain by CMD. This density is equal to nmd/ntmd where nmd, ntmd is the number of mobile atoms and total atoms in the mobile domain, respectively. Because of immobile atoms are concentrated in immobile domains and most of mobile atoms reside in the mobile domain, we have: CIMD > CCS and CMD > CCS. This leads to that large SixOy subnets of immobile and mobile atom are formed in the separate domains as observed from the simulation. Moreover, due to CIMD > CMD the number of bonds for immobile atoms is larger than mobile atoms. With increasing temperature, the size of immobile domains decreases and CMD approaches to CCS. To get more details about spatial distributions of selective atoms we calculate the size distribution of link-clusters with rlk = 1.9 Å. The size of link-cluster is defined as a number of atoms belonging to this cluster. The result is presented in Tables 1 and 2. At 3000 K the size of linkcluster for SRA varies from 1 to 6, whereas SMA and SIMA have numerous clusters with size from 7 to 36. The average size of link-cluster for SIMA is significantly larger than that for SMA. At 3500 K those distributions are similar for SRA and SMA which means that mobile atoms are uniformly distributed in the space at high temperature. Moreover they show a clear difference between SMA and SIMA. However, the difference is less pronounced compared to the case of 3000 K. The result seen in Tables 1 and 2 can be interpreted by that immobile atoms are concentrated in immobile domains which leads to forming large link-clusters of immobile atom. In contrast, the majority of mobile atoms reside in the mobile domain so that large link-clusters of mobile atom are present in the mobile domain. The spatial distribution of selective and random atoms can be described schematically in Fig. 6. It can be seen that there are numerous large link-clusters of selective atoms, while random atoms are distributed in the space by small clusters. The system contains only small clusters of random atoms because of CIMD > CCS and CMD > CCS. The average size of link-cluster of immobile atom is significantly larger than the one of mobile atom due to CIMD > CMD. The size distribution of link-cluster is almost unchanged with time because of CIMD and CMD are slightly fluctuated during the time tobs. At 3500 K the size distribution of link-cluster is similar for SRA and SMA due to that the average size of immobile
Table 1 The number of link-clusters in configurations at t = 71.7 and 143.4 ps, and at 3000 K. Here SCl, NCl is the size of link-cluster and number of clusters having the corresponding SCl, respectively. SCl
1 2 3 4 5 6 7 8 9 11 16 17 19 20 36 –
NCl, t = 71.7 ps
SCl
SMA
SIMA
SRA
310 14 22 10 8 4 3 3 1 2 1 0 0 0 0 –
135 61 18 8 6 5 2 8 3 0 0 1 1 1 1 –
338 75 22 6 2 2 0 0 0 0 0 0 0 0 0 –
1 2 3 4 5 6 7 8 9 10 11 13 14 16 23 30
NCl, t = 143.4 ps SMA
SIMA
SRA
312 27 23 9 2 4 2 2 1 2 0 1 0 0 1 0
149 49 17 10 9 5 2 8 2 2 1 0 1 1 0 1
345 68 26 5 3 1 0 0 0 0 0 0 0 0 0 0
Table 2 The number of link-clusters in the configuration at t = 71.7 and 143.4 ps, and at 3500 K. SCl
1 2 3 4 5 6 7 8
NCl, t = 71.7 ps
SCl
SMA
SIMA
SRA
396 30 24 15 1 0 1 0
268 59 28 16 5 2 3 1
335 76 22 6 2 1 1 0
1 2 3 4 5 6 7 8
NCl, t = 143.4 ps SMA
SIMA
SRA
395 37 20 5 4 4 1 0
307 65 14 9 6 4 1 3
361 61 29 5 2 0 0 0
Fig. 6. The spatial distribution of random and selective atoms. Here we represent the atom by a circle with radius of rlk / 2. If two atoms form a linkage, then the circles overlap with each other. A linkcluster therefore corresponds to a subset of overlapped circles. The red, yellow and blue circle represents the random, immobile and mobile atom, respectively. The other-type atoms are not shown here. It can be seen that selective atoms form numerous large linkclusters, whereas random atoms form only small ones. The average size of link-cluster of immobile atoms is larger than that of mobile atoms. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
4
Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx
P.K. Hung et al.
Table 3 The time evolution of link-cluster with size larger than four for SMA at 3000 K. t=0
t = 14.34 ps
t = 71.7 ps
t = 143.4 ps
SCl
NCl
SCl
NCl
SCl
NCl
SCl
NCl
5 6 7 8 9 12 15 – – – – –
6 1 2 3 2 1 1 – – – – –
1 3 4 5 6 7 8 9 12 16 – –
3 1 3 4 2 3 2 2 1 1 – –
1 2 3 4 5 6 7 10 11 12 – –
15 5 7 7 3 2 3 1 2 2 – –
1 2 3 4 5 6 7 8 9 10 13 23
24 9 8 4 2 3 1 1 1 2 1 1
Table 4 The time evolution of link-cluster with size larger than four for SIMA at 3000 K. t=0
t = 14.34 ps
t = 71.7 ps
Fig. 7. The number of defective atoms as a function of time.
t = 143.4 ps
SCl
NCl
SCl
NCl
SCl
NCl
SCl
NCl
5 6 7 8 9 11 17 18 – – – – – – –
7 7 2 7 3 1 2 2 – – – – – – –
1 2 3 4 5 6 7 8 9 10 11 12 17 18 20
3 1 1 3 6 5 1 7 4 1 1 1 1 1 1
1 2 3 5 6 7 8 9 10 14 17 19 34 – –
3 1 1 6 3 2 9 2 1 1 1 2 1 – –
1 2 3 4 5 6 7 8 9 10 11 14 16 30 –
6 4 1 2 6 5 2 8 2 2 1 1 1 1 –
of the found list and may consist of other selective atoms. Finally, we determine reforming clusters in different configurations at t > 0. The calculation result for four configurations is shown in Tables 3 and 4. We note that there are a number of small reforming clusters with size less than four. For SMA the number of small cluster largely increases with time (see Table 5). This indicates that initial clusters strongly split into small clusters. Therefore the atoms make a strong rearrangement in the mobile domain. On other hands, the number of large clusters with size bigger than four slightly decreases from 16 to 12. The number of total atoms belonging to large clusters also weakly changes. In particular, NTA varies from 104 to 115 that slightly smaller than the one of initial clusters equal to 119. From this follows that initial clusters are split into small clusters which then merge with other atoms to form large clusters with size up to 23 atoms. This evidences that CMD is much larger than CCS. It is expected that the mobile atoms are redistributed from mobile to immobile domains due to CIMD > CMD. Consequently, CMD and the average size of reforming clusters decreases as t increases. However, the average size of reforming clusters NTA/NLCl slightly varies during the time tobs. This evidences that the exchange of mobile atoms between mobile and immobile domains proceeds so that CMD is slightly fluctuated. In the case of SIMA the number of small cluster increases from 8 to 13 that significantly smaller compared to the case of SMA. There are numerous large clusters with size up to 34 atoms. NTA/NLCl is also slightly changes with time. This result evidences CIMD slightly varies during the time tobs. Also the rate of atomic rearrangement in the immobile domain is significantly smaller than the one in the mobile domain. We calculate the number of defective atoms from considered sets. The defective atom is defined as over-coordinated Si or O. Most of defective atoms have the coordination number of five for Si and three for O. The Si with the coordination number of four and O having two coordinated atoms is called normal atom. In Fig. 7 we plot the number of defective atom as a function of time. At 3000 K the number of defective atoms of SMA is significantly larger than the one of SIMA. This is originated from that the diffusion is mainly realized by reactions (1) via defective units such as SiO5 and OSi3. We conclude that the immobile domains have large fraction of normal atom. In contrast, the mobile domain possesses large fraction of defective atoms. Therefore the immobile and mobile domains possess different local microstructures
domains significantly decreases and CMD is close to CCS. This result indicates that the liquid exhibits DH the degree of which reduces with increasing temperature. To get more information about mobile and immobile domains we examine the evolution of large link-clusters with rlk = 1.9 Å. The calculation is performed as follows. Firstly, we determine all link-clusters with size larger than four in the starting configuration (t = 0). The clusters found are called initial clusters. Next, a list of atoms belonging to initial clusters is determined. As the time proceeds, the atoms of the found list are redistributed through new link-clusters called reforming clusters. Unlike initial cluster, the reforming cluster has at least an atom Table 5 Characteristics of initial and reforming clusters. Here NSCl, NLCl is the number of cluster with size SCl < 5 and SCl ≥ 5, respectively; NTA is the number of total atoms belonging to the clusters with size SCl ≥ 5. t, ps
0 14.34 71.7 143.4
SMA
SIMA
NSCl
NLCl
NTA
NSCl
NLCl
NTA
0 7 34 45
16 15 13 12
119 115 104 108
0 8 11 13
31 29 28 29
255 247 265 247
5
Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx
P.K. Hung et al.
Fig. 8. The link-cluster function of SMA and SIMA for tobs = 71.7 and 143.4 ps at 3000 (a) and 3500 K (b).
two-domain structure is observed within a time interval [0, tobs]. For the next interval [tobs, 2tobs] we observe other immobile domains which possess similar properties. In particular, the size distribution of immobile domain is almost unchanged, but their shape and position change significantly. Overall, the scenario for dynamics can be described as follows. During a moderate time tobs the immobile atoms are concentrated in immobile domains, which lead to forming large link-clusters with rlk = 1.9 Å. Meanwhile most of mobile atoms move inside a mobile domain so that large link-clusters of mobile atom are formed there. The rate of atomic rearrangement in the immobile domain is significantly smaller than the one in the mobile domain. The exchange of selective
which mean that DH is accompanied with the structure heterogeneity in the liquid. In the case of 3500 K the number of defective atoms is close for SIMA and SMA. This result can be interpreted by that at this temperature the rate of reaction (1) is high and the spatial distribution of selective atoms is more homogeneous. In Fig. 8 we present Flink(r, t) for different tobs. One can see that Flink(r, t) for selective atoms approaches to the one for random atoms as the time tobs increases. This means that the spatial correlation between selective atoms is weaker as tobs is longer. As a result, the size of immobile domains decreases with increasing tobs, i.e. upon long time tobs we do not observe the two-domain structure. It is worth to note that the
6
Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx
P.K. Hung et al.
the spatial distribution of selective atoms is non-uniform, the reactions are also non-uniformly distributed in the space. Due to the SieO bond the atoms do not make random displacements, but they tend to move by a group as a like-molecule SixOy, where x > 0, y > 0. For instance, Si2O7 kept eight init-bonds. The SixOy molecule with lifetime t is defined as a subset of atoms connected to each other via init-bonds over whole time t. Table 6 presents the number of such molecules for different lifetimes. As expected, the size of large SixOy decreases as the time t increases. We observed that there are numerous large SixOy molecules existed for long time even at high temperature. This clearly demonstrates the cooperative movement of atoms via SixOy. Such movement results in that the selective atoms tend to locate nearby to form Si-O bonds. Moreover the selective atoms form numerous large link-clusters with rlk = 1.9 Å. We conclude that DH observed is accompanied with the cooperative movement via SixOy. 4. Conclusions MD simulation is carried out to study the diffusion mechanism, dynamics and structure heterogeneity in the liquid silica. We propose link-cluster and init-bond functions to explore these phenomena. The calculation shows that at low temperature selective atoms are nonuniformly distributed in the space. In particular, the selective atoms form numerous large link-clusters with rlk = 1.9 Å. At 3000 K the immobile and mobile atoms form respectively large clusters with size up to 34 and 23 atoms. Whereas the random atoms are arranged by small clusters with size < 6 atoms. We reveal that during a moderately long time tobs the liquid has a two-domain structure that consists of separate immobile and mobile domains. Immobile domains strongly differ from the mobile domain in the rate of atomic rearrangement and that the fraction of defective atoms in immobile domains is significantly larger than the one in the mobile domain. The calculation also reveals that the diffusion is realized via bondbreaking mechanism. The reactions are non-uniformly distributed in the space and the liquid exhibits DH the degree of which reduces with increasing temperature. Furthermore the simulation shows a strong correlation between mobility of atoms and init-bond function. We find that DH is accompanied with structure heterogeneity and cooperative movement of atoms via like-molecules SixOy.
Fig. 9. The Finb(r) as a function of time for SMA, SIMA and SRA.
Table 6 The number of SixOy like-molecules. Smolecule, Nmolecule is the range of molecule size and number of molecules having the corresponding sizes, respectively; T is the temperature. t = 71.7 ps; T = 3000 K
t = 143.4 ps; T = 3000 K
t = 71.7 ps; T = 3500 K
t = 143.4 ps; T = 3500 K
Smolecule
Nmolecule
Smolecule
Nmolecule
Smolecule
Nmolecule
Smolecule
Nmolecule
2–9 10–12 2469
91 3 1
2–9 10–37 1549
248 11 1
2–9 10–14 20
483 3 1
2 3 4–5
149 25 4
Acknowledgement
atoms between mobile and immobile domains proceeds so that CIMD and CMD are slightly fluctuated within the time tobs. We calculate Finb(t) by averaging over all O and Si which belong to SIMA, SMA and SRA. The result is plotted in Fig. 9. Clearly that the mobility of atom is strongly correlated with Finb(t). Namely, Finb(t) decreases in the order: SIMA → SMA → SRA. Also at 3500 K the decrease of Finb(t) with time proceeds much faster than at 3000 K. These results support the bond-breaking mechanism by the fact that the atoms possessing high mobility undergo more reactions, i.e. more init-bonds are replaced. The fast decrease of Finb(t) with time at 3500 K is caused by that the SieO bond is broken more easy at high temperature. Because
The authors are grateful for support by the NAFOSTED Vietnam (grant 103.01-2015.12). This research was supported in part by MEXT as “Exploratory Challenge on Post-K computer” (Frontiers of Basic Science: Challenging the Limits) and by RIKEN as the iTHES Project. References [1] H.C. Andersen, Proc. Natl. Acad. Sci. U. S. A. 102 (2005) 6686. [2] F. Sciortino, P. Tartaglia, Phys. 54 (2005) 471.
7
Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx
P.K. Hung et al.
[3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
[15] Sharon C. Glotzer, J. Non-Cryst. Solids 274 (2000) 342. [16] T. Kawasaki, T. Araki, H. Tanaka, Phys. Rev. Lett. 99 (2007) 215701. [17] Wen-Sheng Xu, Zhao-Yan Sun, Li-Jia An, J. Phys. Condens. Matter 24 (2012) 325101. [18] G.N. Greaves, S. Sen, Adv. Phys. 56 (2007) 1. [19] J.S. Langer, Rep. Prog. Phys. 77 (2014) 042501. [20] T. Kawasaki, H. Tanaka, J. Phys. Condens. Matter 23 (2011) 194121. [21] E.R. Weeks, J.C. Crocker, D.A. Weitz, J. Phys. Condens. Matter 19 (2007) 205131. [22] E.R. Weeks, D.A. Weitz, Phys. Rev. Lett. 89 (2002) 095704. [23] C. Donati, S.C. Glotzer, P.H. Poole, Phys. Rev. Lett. 82 (1999) 5064. [24] D. Coslovich, G. Pastore, J. Phys. Condens. Matter 21 (2009) 285107. [25] Y. Wang, et al., Nat. Commun. 5 (2014) 3241. [26] B. van Beest, G. Kramer, R. van Santen, Phys. Rev. Lett. 64 (1990) 1955. [27] P.K. Hung, L.T. Vinh, J. Non-Cryst. Solids 352 (2006) 5531. [28] P.K. Hung, N.T.T. Ha, N.V. Hong, Eur. Phys. J. E 36 (2013) 60.
R. Richert, J. Phys. Condens. Matter 14 (2002) 703. T. Kawasaki, H. Tanaka, J. Phys. Condens. Matter 22 (2010) 232102. I. Chang, H. Sillescu, J. Phys. Chem. B 101 (1997) 8794. H. Sillescu, J. Non-Cryst. Solids 243 (1999) 81. M.T. Cicerone, M.D. Ediger, J. Chem. Phys. 104 (1996) 7210. B. Schiener, R.V. Chamberlin, G. Diezemann, R. Bohmer, J. Chem. Phys. 107 (1997) 7746. R. Richert, J. Phys. Chem. B 101 (1997) 6323. M. Shajahan, G. Razul, Gurpreet S. Matharoo, Peter H. Poole, J. Phys. Condens. Matter 23 (2011) 235103. Richard K. Darst, David R. Reichman, Giulio Biroli, J. Chem. Phys. 132 (2010) 044510. Z. Zheng, et al., Nat. Commun. 5 (2014) 3829. G.A. Appignanesi, J.A. Rodriguez Fris, J. Phys. Condens. Matter 21 (2009) 203103. Tanaka, et al., Nat. Mater. 9 (2010) 324.
8