Properties of a continuous-random-network model for amorphous systems

Yuhai Tu, J. Tersoff, G. Grinstein IBM Research Division, T. J. Watson Research Center, P. O. Box 218, Yorktown Heights, NY 10598    David Vanderbilt Department of Physics, Rutgers University, P. O. Box 849, Piscataway NJ 08855
July 16, 2016
Abstract

We use a Monte Carlo bond-switching method to study systematically the thermodynamic properties of a “continuous random network” model, the canonical model for such amorphous systems as a-Si and a-SiO 22{}_{2} . Simulations show first-order “melting” into an amorphous state, and clear evidence for a glass transition in the supercooled liquid. The random-network model is also extended to study heterogeneous structures, such as the interface between amorphous and crystalline Si.

pacs:
PACS numbers: 61.43.-j, 64.70.Pf

Amorphous materials have been intensely studied because of both their technological importance and fundamental interest. Due to the complexity of these systems, theoretical studies have relied mainly on numerical simulations. The standard method is molecular dynamics (MD) (1); (2); (3); (4), where atoms are represented by point particles with more or less realistic interactions, and the equations of motion are integrated numerically.

Most studies of glassy behavior in amorphous systems have involved hard-sphere or Lennard-Jones models for metallic glasses. However, the most important and ubiquitous amorphous materials are network glasses such as a-SiO 22{}_{2} . These may also be studied by MD, but large-scale or long-time simulations are limited to the few systems where suitable classical models are available for the atomic interactions (1); (4). Moreover, MD methods are inherently less efficient for strong network-forming materials, because there is a large energy barrier to breaking and reforming bonds in these systems. Most of the computer time is therefore spent in following local vibrations, waiting for the infrequent bond-switching events.

An alternative is to use the canonical model for such network glasses: a continuous random network (CRN) of atoms and bonds (5). The most realistic available models of amorphous Si are of this type. They have been generated by Wooten, Winer and Weaire (6)(hereafter WWW), using an ingenious Monte Carlo (MC) approach for generating random networks by bond switching.

Here we show that this approach can be extended to directly study the thermodynamic properties of disordered materials. The system is represented explicitly as a CRN, i.e., as a set of atoms and a “neighbor list” specifying which pairs of atoms are connected by bonds. Equilibrium or quasi-equilibrium properties are determined by MC sampling, using the bond-switching move introduced by WWW.

The specific model system studied here has four bonds to each atom, as for a-Si. But we view it principally as a generic CNR model, and we systematically study its phase diagram and glass transition. This approach, and the phase structure we find, should apply also to more complicated covalent amorphous solids such as a-SiO 22{}_{2} . In addition, we use the model to explore the properties of the interface between the crystalline and amorphous phases of silicon.

The energy of the CRN is given by a Keating-like valence force model (7), and depends on both the positions {ri}subscript𝑟𝑖\{\vec{r}_{i}\} of the atoms and the set ψ𝜓\psi of bonds connecting pairs of atoms:

Etot(ψ,{ri})subscript𝐸tot𝜓subscript𝑟𝑖\displaystyle E_{\rm tot}(\psi,\{\vec{r}_{i}\}) =\displaystyle= i,jϵψ12kθb02(cosθij-cosθ0)2subscript𝑖𝑗italic-ϵ𝜓12subscript𝑘𝜃superscriptsubscript𝑏02superscriptsubscript𝜃𝑖𝑗subscript𝜃02\displaystyle\sum_{i,j\epsilon\psi}\frac{1}{2}k_{\theta}b_{0}^{2}(\cos\theta_{% ij}-\cos\theta_{0})^{2} (1)
+jϵψ12kb(bj-b0)2.subscript𝑗italic-ϵ𝜓12subscript𝑘𝑏superscriptsubscript𝑏𝑗subscript𝑏02\displaystyle+\sum_{j\epsilon\psi}\frac{1}{2}k_{b}(b_{j}-b_{0})^{2}~{}.

Here j𝑗j represents the j𝑗j th bond, bjsubscript𝑏𝑗b_{j} is its length, θijsubscript𝜃𝑖𝑗\theta_{ij} is the angle between bonds i𝑖i and j𝑗j connected to a common atom, b0subscript𝑏0b_{0} is the preferred bond length, θ0subscript𝜃0\theta_{0} is the preferred bond angle, and kθsubscript𝑘𝜃k_{\theta} and kbsubscript𝑘𝑏k_{b} are “spring constants.In order to focus on the role of network structure, we write the energy as a function solely of bond topology, minimizing Etotsubscript𝐸totE_{\rm tot} with respect to the geometrical coordinates {ri}subscript𝑟𝑖\{\vec{r}_{i}\} :

E(ψ)=min{ri}Etot(ψ,{ri}).𝐸𝜓subscriptsubscript𝑟𝑖subscript𝐸tot𝜓subscript𝑟𝑖E(\psi)=\min_{\{\vec{r}_{i}\}}E_{\rm tot}(\psi,\{\vec{r}_{i}\})~{}~{}. (2)

Using E(ψ)𝐸𝜓E(\psi) , we can study the statistical properties of the system through MC simulation. We use the WWW construction for the local MC moves in the ψ𝜓\psi space. That is, from an initial configuration ψ1subscript𝜓1\psi_{1} , a bond is chosen randomly (call it BC), and one more bond connected to each terminus is also chosen randomly (bonds BA and CD), the only constraint being that all four atoms A, B, C and D must be distinct. The switching move leading to ψ2subscript𝜓2\psi_{2} is then simply the cutting of bonds BA and CD and the formation of new bonds AC and BD. In this way the system samples topologically distinct configurations without introducing “dangling bonds” or changing the number of bonds to any atom.

For a given temperature T𝑇T , one uses the Metropolis MC algorithm to decide whether a given switching move is accepted or rejected: the move is accepted with probability P=min[1,exp(-ΔE/kBT)]𝑃1Δ𝐸subscript𝑘𝐵𝑇P=\min[1,\exp(-\Delta E/k_{B}T)] , where ΔE=E(ψ2)-E(ψ1)Δ𝐸𝐸subscript𝜓2𝐸subscript𝜓1\Delta E=E(\psi_{2})-E(\psi_{1}) is the energy change. Note that the temperature T𝑇T applies only to the bond topology; the system is always in its ground state with respect to phonons. This makes the model less realistic in describing certain aspects of specific systems, but more generic in distilling the role of network topology and excluding other issues.

We choose parameter values appropriate for a-Si: kθ=0.647subscript𝑘𝜃0.647k_{\theta}=0.647 eV/Å 22{}^{2} , kb=9.08subscript𝑘𝑏9.08k_{b}=9.08 eV/Å 22{}^{2} , b0=2.35subscript𝑏02.35b_{0}=2.35 Å, and θ0=109subscript𝜃0superscript109\theta_{0}=109^{\circ} . In addition, it is important that the neighbor list ψ𝜓\psi remain consistent with the geometry {ri}subscript𝑟𝑖\{\vec{r}_{i}\} , with close atoms included in the neighbor list, and distant atoms excluded. To guarantee this, we include an extra energy term Esuperscript𝐸E^{\prime} in Etotsubscript𝐸totE_{\rm tot} to prevent structures with “false” neighbors from occurring: E=γmn(d2-|rm-rn|)3superscript𝐸𝛾subscript𝑚𝑛superscriptsubscript𝑑2subscript𝑟𝑚subscript𝑟𝑛3E^{\prime}=\gamma\sum_{mn}(d_{2}-|\vec{r}_{m}-\vec{r}_{n}|)^{3} . Here m𝑚m and n𝑛n label atoms which are neither 1st nor 2nd neighbors in ψ𝜓\psi , but for which |rm-rn|subscript𝑟𝑚subscript𝑟𝑛|\vec{r}_{m}-\vec{r}_{n}| is actually less than the distance d2subscript𝑑2d_{2} =3.84 Å between next-nearest neighbors in crystalline silicon. We use γ=0.5𝛾0.5\gamma=0.5 eV/Å 33{}^{3} .

We begin our simulations in the ground state of the system, a diamond-structure crystal, where each atom has four bonds. At low temperature this crystalline phase is stable, while at high temperature the crystal “melts” into a disordered “liquid” phase. Note that this disordered system is still a network, with the same bond coordination as the crystal. Good glass-formers such as SiO 22{}_{2} typically share this characteristic that the liquid retains the network structure (3).

It is difficult to determine the melting point accurately in a homogenous system, because of the energy barrier to nucleation of a new phase. We therefore create a system with a solid-liquid interface, and study the interface motion as a function of temperature. For temperatures T<Tm𝑇subscript𝑇𝑚T<T_{m} , the crystal phase will invade the amorphous phase, while for T>Tm𝑇subscript𝑇𝑚T>T_{m} the reverse occurs. The transition temperature Tmsubscript𝑇𝑚T_{m} corresponds to the temperature where the interface motion vanishes.

The sample is prepared by taking a crystal and allowing WWW switching in half the cell, while forbidding switching in the other half. At sufficiently high T𝑇T , the first half becomes liquid. To determine the interface motion, the system is then cooled to the temperature of interest, and switching is allowed throughout the sample.

We have simulated interface motion in this manner at different temperatures for systems with N𝑁N =432 atoms, with interfaces oriented along the [111] direction. Figure 1 shows snapshots of the atomic positions projected onto the (011¯)01¯1(01\bar{1}) plane, at six successive times, for each of the two temperatures kBT1subscript𝑘𝐵subscript𝑇1k_{B}T_{1} =0.50 eV [Fig. 1(a)] and kBT2subscript𝑘𝐵subscript𝑇2k_{B}T_{2} =0.55 eV [Fig. 1(b)]. Fig. 1 makes clear that the amorphous central region shrinks at temperature T1subscript𝑇1T_{1} by recrystallization, while at T2subscript𝑇2T_{2} the more stable amorphous phase grows at the expense of the crystal. We conclude that T1<Tm<T2subscript𝑇1subscript𝑇𝑚subscript𝑇2T_{1}<T_{m}<T_{2} . From further calculations of this type, we estimate that kBTm0.53±0.01similar-tosubscript𝑘𝐵subscript𝑇𝑚plus-or-minus0.530.01k_{B}T_{m}\sim 0.53\pm 0.01 eV.

This Tmsubscript𝑇𝑚T_{m} is a very high temperature, around 6000K. Note that real Si does not form such a network liquid, but rather melts into a very different high-coordination liquid phase at much lower temperature, around 1700K. As a result, Si is a poor glass-former, and a-Si can only be formed by processes which are very far from equilibrium.

We have also studied the properties of the homogeneous disordered network (the “liquid” phase) for temperatures both above and below Tmsubscript𝑇𝑚T_{m} , for system size N𝑁N =216. In particular, the average energy per atom, Ea(T)subscript𝐸𝑎𝑇E_{a}(T) , of the amorphous phase is determined by equilibrating for 2×105similar-toabsent2superscript105\sim 2\times 10^{5} MC steps and then averaging over 106similar-toabsentsuperscript106\sim 10^{6} MC steps, for each T𝑇T . To further reduce inaccuracy due to statistical fluctuation, the result is averaged over typically 10-20 runs with different random number seeds. From standard thermodynamic relations, we can then easily calculate the entropy, Sa(T)subscript𝑆𝑎𝑇S_{a}(T) , and free energy, Fa(T)subscript𝐹𝑎𝑇F_{a}(T) , per atom: Sa(T)=Sa(T0)+T0T(1/T)(Ea/T)𝑑Tsubscript𝑆𝑎𝑇subscript𝑆𝑎subscript𝑇0superscriptsubscriptsubscript𝑇0𝑇1𝑇subscript𝐸𝑎𝑇differential-d𝑇S_{a}(T)=S_{a}(T_{0})+\int_{T_{0}}^{T}(1/T)(\partial E_{a}/\partial T)dT , and Fa(T)=Ea(T)-TSa(T)subscript𝐹𝑎𝑇subscript𝐸𝑎𝑇𝑇subscript𝑆𝑎𝑇F_{a}(T)=E_{a}(T)-TS_{a}(T) . It is particularly convenient to choose the arbitrary temperature T0subscript𝑇0T_{0} to be Tmsubscript𝑇𝑚T_{m} , since Fa(Tm)=Fc(Tm)subscript𝐹𝑎subscript𝑇𝑚subscript𝐹𝑐subscript𝑇𝑚F_{a}(T_{m})=F_{c}(T_{m}) , where Fc(Tm)subscript𝐹𝑐subscript𝑇𝑚F_{c}(T_{m}) is the free energy for the crystal. In this model Fc(T)=0subscript𝐹𝑐𝑇0F_{c}(T)=0 to an excellent approximation for temperatures in the range of interest (8).

The resulting values for Ea(T)subscript𝐸𝑎𝑇E_{a}(T) , Fa(T)subscript𝐹𝑎𝑇F_{a}(T) , and Sa(T)subscript𝑆𝑎𝑇S_{a}(T) are shown in Fig. 2. When T<Tm𝑇subscript𝑇𝑚T<T_{m} , Fa(T)>0subscript𝐹𝑎𝑇0F_{a}(T)>0 , so the crystal phase is thermodynamically preferred, while for T>Tm𝑇subscript𝑇𝑚T>T_{m} , Fa(T)<0subscript𝐹𝑎𝑇0F_{a}(T)<0 , and the amorphous “liquid” phase is more stable. Our simulations clearly indicate that the crystalline phase is metastable in the region of stability of the amorphous phase, and vice versa, with a large nucleation barrier separating the phases.

Thus, we are able to supercool the liquid and obtain well-defined quasi-equilibrium properties for the metastable liquid below the first order transition at Tmsubscript𝑇𝑚T_{m} . The energy and entropy curves for this liquid exhibit fairly abrupt reproducible changes in slope at a rather well-defined temperature, kBTg0.4similar-tosubscript𝑘𝐵subscript𝑇𝑔0.4k_{B}T_{g}\sim 0.4 eV. This suggests that the liquid phase undergoes a glass transition at kBTg0.4similar-tosubscript𝑘𝐵subscript𝑇𝑔0.4k_{B}T_{g}\sim 0.4 eV, so Tg0.75Tmsimilar-tosubscript𝑇𝑔0.75subscript𝑇𝑚T_{g}\sim 0.75T_{m} .

It is interesting to note in Fig. 2 that extrapolations of Ea(T)subscript𝐸𝑎𝑇E_{a}(T) and Sa(T)subscript𝑆𝑎𝑇S_{a}(T) from above Tgsubscript𝑇𝑔T_{g} would give negative values at very low temperature. This phenomenon is known as the Kauzmann paradox, and is rather common for materials that form structural glasses (9). Of course Fig. 2 shows that the energy and entropy for kBT<0.4eVsubscript𝑘𝐵𝑇0.4𝑒𝑉k_{B}T<0.4eV lie significantly above the curves extrapolated from higher T𝑇T , thereby avoiding a true, unphysical paradox.

One can study this transition more quantitatively by calculating the “time”-dependent structure factor (10)

S(q,t)=ρ(q,t)ρ(-q,t+t),𝑆𝑞𝑡delimited-⟨⟩𝜌𝑞superscript𝑡𝜌𝑞𝑡superscript𝑡S(\vec{q},t)=\langle\rho(\vec{q},t^{\prime})\rho(-\vec{q},t+t^{\prime})\rangle% ~{}~{}, (3)

where ρ(q,t)𝜌𝑞superscript𝑡\rho(\vec{q},t^{\prime}) is the Fourier transform of the atomic density at time tsuperscript𝑡t^{\prime} , ρ(q,t)=j=1Nexp(iqrj(t))/N𝜌𝑞superscript𝑡superscriptsubscript𝑗1𝑁𝑖𝑞subscript𝑟𝑗superscript𝑡𝑁\rho(\vec{q},t^{\prime})=\sum_{j=1}^{N}\exp(i\vec{q}\cdot\vec{r}_{j}(t^{\prime% }))/\sqrt{N} , and the angle brackets denote an average over tsuperscript𝑡t^{\prime} . In the liquid phase, S(q,t)𝑆𝑞𝑡S(\vec{q},t) decays exponentially, S(q,t)exp(-D|q|2t)similar-to𝑆𝑞𝑡𝐷superscript𝑞2𝑡S(\vec{q},t)\sim\exp(-D|\vec{q}|^{2}t) , where D𝐷D is the diffusion constant of the liquid (11). In a glass phase, however, where there is frozen (albeit random) structure, S(q,t)𝑆𝑞𝑡S(\vec{q},t) in the thermodynamic limit remains finite as t𝑡t\rightarrow\infty , decaying only because of the finite size of the simulated system. We choose |q|=2π/b0𝑞2𝜋subscript𝑏0|\vec{q}|=2\pi/b_{0} , with b0subscript𝑏0b_{0} the bond length of the crystal structure, so that S(q,t)𝑆𝑞𝑡S(\vec{q},t) will decay as rapidly as possible in the liquid phase. In addition to averaging over tsuperscript𝑡t^{\prime} , we average over many orientations of q𝑞\vec{q} to get better statistics. The results are shown in Fig. 3, where the normalized time-dependent structure factor S(|q|,t)/S(|q|,0)𝑆𝑞𝑡𝑆𝑞0S(|\vec{q}|,t)/S(|\vec{q}|,0) is plotted versus time delay t𝑡t , for four different temperatures around Tgsubscript𝑇𝑔T_{g} . (The “time” units here are MC steps per atom.)

It can be seen from Fig. 3 that the structure factor decays to zero for kBTsubscript𝑘𝐵𝑇k_{B}T =0.40 and 0.45 eV, but remains finite for the lower temperatures, kBTsubscript𝑘𝐵𝑇k_{B}T =0.30 and 0.35 eV. Of course the number of accepted MC switching steps decreases with decreasing temperature. However, plotting the structure factor vs. accepted steps gives a similar picture, leading to the same conclusion. The simulations ran long enough to produce more than 2000 accepted switching moves even for the lowest temperature, kBTsubscript𝑘𝐵𝑇k_{B}T =0.3 eV.

These results confirm the conclusion above, that there is a well-defined glass transition in the idealized CRN model (12), and this transition occurs at kBTg0.40similar-tosubscript𝑘𝐵subscript𝑇𝑔0.40k_{B}T_{g}\sim 0.40 eV for the specific network considered here. Thus the CRN model has the following phase structure: In equilibrium, as temperature decreases there is a first-order transition from a disordered network “liquid” phase to the crystalline phase. The liquid phase may be supercooled below Tmsubscript𝑇𝑚T_{m} , and at T=Tg<Tm𝑇subscript𝑇𝑔subscript𝑇𝑚T=T_{g}<T_{m} undergoes a transition into a metastable glassy phase.

It is interesting to note that the model glass phase obtained by cooling the metastable liquid through Tgsubscript𝑇𝑔T_{g} has essentially the same structure as was shown by WWW (6)to provide the best available model for the structure of real a-Si. Yet real a-Si is formed by radically different processes, very far from equilibrium. This suggests the existence of a rather well-defined metastable amorphous phase, whose structure (after annealing) is relatively independent of the kinetic history of the material.

As demonstrated in our calculation of Tmsubscript𝑇𝑚T_{m} , the random network model can easily be adapted to study the interface between amorphous and crystalline phases. To investigate the effects of crystalline anisotropy, we prepared interfaces as before, oriented parallel to the (100), (110), and (111) planes of the crystal. We then equilibrated at T=Tm𝑇subscript𝑇𝑚T=T_{m} , the only temperature for which the interface is stationary. After initial transients, the (100) and (110) interfaces develop {111} facets, while the (111) interface remains flat. For the (100) interface, Fig. 4(a) shows the positions of all the Si atoms at one instant of time. It is clear that the interface is unstable the dotted lines indicate two {111} facets which have formed. Similar faceting occurs for the (110) interface. Fig. 4(b) shows that, in contrast, the (111) interface remains planar. Whether such faceting occurs in real a-Si will depend on the relative importance of thermodynamic and kinetic factors.

We can study the decay of crystalline order near the interface in the amorphous phase by calculating the local energy density, which is zero for the perfect crystal and non-zero in the amorphous phase. We find that the partial crystalline structure on the amorphous side of the interface always forms in double layers, consistent with the crystal structure in the [111] direction. This partial order decays rapidly as one moves away from the interface, with a decay length of roughly one (111) double layer.

In summary, we have studied the properties of a random-network model for amorphous materials. While we have used parameters for a-Si, the general behavior should be rather generic. We demonstrated the existence of a glass transition in an ideal network model. We have also applied the model to study the crystal-amorphous interface. For Si the (111) interface was found to be stable, but (110) and (001) interfaces are unstable against formation of (111) facets. This general approach should be applicable to a variety of systems which maintain a nearly ideal network structure, as well as for studies of generic aspects of network glasses.

D.V. acknowledges support of NSF Grant DMR-9613648.

References

  • (1) K. Ding and H. C. Andersen, Phys. Rev. B 34 , 6987(1986); R. Biswas, G. S. Grest and C. M. Soukoulis, Phys. Rev. B 36 , 7437 (1987); W. D. Luedtke and U. Landman, Phys. Rev. B 37 , 4656 (1988).
  • (2) R. Car and M. Parrinello, Phys. Rev. Lett. 60 , 204 (1988).
  • (3) J. Sarnthein, A. Pasquarello, and R. Car, Phys. Rev. B 52 , 12690 (1995).
  • (4) K. Vollmayr, W. Kob and K. Binder, Phys. Rev. B 54 , 15809 (1996).
  • (5) W. H. Zachariasen, J. Ann. Chem. Soc. 54 , 3841 (1932).
  • (6) F. Wooten, K. Winer, and D. Weaire, Phys. Rev. Lett. 54 , 1392 (1985).
  • (7) P.N. Keating, Phys. Rev. 145 , 637(1966).
  • (8) Since the phonon temperature is always zero, deviations from Fc(T)=0subscript𝐹𝑐𝑇0F_{c}(T)=0 occur only because of an exponentially small number of thermal lattice defects in the crystal.
  • (9) W. Kauzmann, Chem. Rev. 43 , 219 (1948).
  • (10) J.J. Ullo and S. Yip, Phys. Rev. Lett. 54 , 1509(1985).
  • (11) E.g., D. Forster, Hydrodynamic Fluctuations, Broken Symmetry and Correlation Functions (Benjamin, New York, 1975).
  • (12) Using MD, Vollmayr et al. (Ref. (4)) have identified and studied a glass transition in amorphous silica, which also has a random network structure. Consistent with experiment, they find that the glass transition temperature decreases with cooling rate, a phenomenon which we have not investigated, but which likely occurs in our model too. See also R.J. Speedy and P.G. Debenedetti, Mol. Phys. 88 , 1293 (1996), and references therein.
Figure 1: Snapshots of projected atom positions at time intervals of 200 MC steps per atom. (a) T=0.5𝑇0.5T=0.5 eV /kBabsentsubscript𝑘𝐵/k_{B} (T<Tm)𝑇subscript𝑇𝑚(T<T_{m}) , where the crystalline phase eventually takes over the whole space. (b) T=0.55𝑇0.55T=0.55 eV /kBabsentsubscript𝑘𝐵/k_{B} (T>Tm)𝑇subscript𝑇𝑚(T>T_{m}) , where the whole system becomes amorphized.
Figure 2: Energy, entropy, and free energy per atom versus temperature. The break in slope for energy and entropy at T=Tg0.4𝑇subscript𝑇𝑔similar-to0.4T=T_{g}\sim 0.4 eV /kBabsentsubscript𝑘𝐵/k_{B} corresponds to a glass transition. ( T01subscript𝑇01T_{0}\equiv 1 eV/kB𝑒𝑉subscript𝑘𝐵eV/k_{B} is introduced so that entropy can be plotted in units of energy.)
Figure 3: Normalized dynamical structure factor versus time (in MC steps per atom) for temperatures around the glass transition.
Figure 4: Interface morphology for different orientations. The (100) interface develops (111) facets (dotted lines). The atom positions are projected onto the (011¯)01¯1(01\bar{1}) plane.