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

###### 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 ${}_{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.PfAmorphous 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 ${}_{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 ${}_{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 $\{\vec{r}_{i}\}$ of the atoms and the set $\psi$ of bonds connecting pairs of atoms:

$\displaystyle E_{\rm tot}(\psi,\{\vec{r}_{i}\})$ | $\displaystyle=$ | $\displaystyle\sum_{i,j\epsilon\psi}\frac{1}{2}k_{\theta}b_{0}^{2}(\cos\theta_{% ij}-\cos\theta_{0})^{2}$ | (1) | ||

$\displaystyle+\sum_{j\epsilon\psi}\frac{1}{2}k_{b}(b_{j}-b_{0})^{2}~{}.$ |

Here $j$ represents the $j$ th bond, $b_{j}$ is its length, $\theta_{ij}$ is the angle between bonds $i$ and $j$ connected to a common atom, $b_{0}$ is the preferred bond length, $\theta_{0}$ is the preferred bond angle, and $k_{\theta}$ and $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 $E_{\rm tot}$ with respect to the geometrical coordinates $\{\vec{r}_{i}\}$ :

$E(\psi)=\min_{\{\vec{r}_{i}\}}E_{\rm tot}(\psi,\{\vec{r}_{i}\})~{}~{}.$ | (2) |

Using $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 $\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 $\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$ , 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(-\Delta E/k_{B}T)]$ , where $\Delta E=E(\psi_{2})-E(\psi_{1})$ is the energy change. Note that the temperature $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_{\theta}=0.647$ eV/Å ${}^{2}$ , $k_{b}=9.08$ eV/Å ${}^{2}$ , $b_{0}=2.35$ Å, and $\theta_{0}=109^{\circ}$ . In addition, it is important that the neighbor list $\psi$ remain consistent with the geometry $\{\vec{r}_{i}\}$ , with close atoms included in the neighbor list, and distant atoms excluded. To guarantee this, we include an extra energy term $E^{\prime}$ in $E_{\rm tot}$ to prevent structures with “false” neighbors from occurring: $E^{\prime}=\gamma\sum_{mn}(d_{2}-|\vec{r}_{m}-\vec{r}_{n}|)^{3}$ . Here $m$ and $n$ label atoms which are neither 1st nor 2nd neighbors in $\psi$ , but for which $|\vec{r}_{m}-\vec{r}_{n}|$ is actually less than the distance $d_{2}$ =3.84 Å between next-nearest neighbors in crystalline silicon. We use $\gamma=0.5$ eV/Å ${}^{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 ${}_{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<T_{m}$ , the crystal phase will invade the amorphous phase, while for $T>T_{m}$ the reverse occurs. The transition temperature $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$ , 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$ =432 atoms, with interfaces oriented along the [111] direction. Figure 1 shows snapshots of the atomic positions projected onto the $(01\bar{1})$ plane, at six successive times, for each of the two temperatures $k_{B}T_{1}$ =0.50 eV [Fig. 1(a)] and $k_{B}T_{2}$ =0.55 eV [Fig. 1(b)]. Fig. 1 makes clear that the amorphous central region shrinks at temperature $T_{1}$ by recrystallization, while at $T_{2}$ the more stable amorphous phase grows at the expense of the crystal. We conclude that $T_{1}<T_{m}<T_{2}$ . From further calculations of this type, we estimate that $k_{B}T_{m}\sim 0.53\pm 0.01$ eV.

This $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 $T_{m}$ , for system size $N$ =216. In particular, the average energy per atom, $E_{a}(T)$ , of the amorphous phase is determined by equilibrating for $\sim 2\times 10^{5}$ MC steps and then averaging over $\sim 10^{6}$ MC steps, for each $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, $S_{a}(T)$ , and free energy, $F_{a}(T)$ , per atom: $S_{a}(T)=S_{a}(T_{0})+\int_{T_{0}}^{T}(1/T)(\partial E_{a}/\partial T)dT$ , and $F_{a}(T)=E_{a}(T)-TS_{a}(T)$ . It is particularly convenient to choose the arbitrary temperature $T_{0}$ to be $T_{m}$ , since $F_{a}(T_{m})=F_{c}(T_{m})$ , where $F_{c}(T_{m})$ is the free energy for the crystal. In this model $F_{c}(T)=0$ to an excellent approximation for temperatures in the range of interest (8).

The resulting values for $E_{a}(T)$ , $F_{a}(T)$ , and $S_{a}(T)$ are shown in Fig. 2. When $T<T_{m}$ , $F_{a}(T)>0$ , so the crystal phase is thermodynamically preferred, while for $T>T_{m}$ , $F_{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 $T_{m}$ . The energy and entropy curves for this liquid exhibit fairly abrupt reproducible changes in slope at a rather well-defined temperature, $k_{B}T_{g}\sim 0.4$ eV. This suggests that the liquid phase undergoes a glass transition at $k_{B}T_{g}\sim 0.4$ eV, so $T_{g}\sim 0.75T_{m}$ .

It is interesting to note in Fig. 2 that extrapolations of $E_{a}(T)$ and $S_{a}(T)$ from above $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 $k_{B}T<0.4eV$ lie significantly above the curves extrapolated from higher $T$ , thereby avoiding a true, unphysical paradox.

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

$S(\vec{q},t)=\langle\rho(\vec{q},t^{\prime})\rho(-\vec{q},t+t^{\prime})\rangle% ~{}~{},$ | (3) |

where $\rho(\vec{q},t^{\prime})$ is the Fourier transform of the atomic density at time $t^{\prime}$ , $\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 $t^{\prime}$ . In the liquid phase, $S(\vec{q},t)$ decays exponentially, $S(\vec{q},t)\sim\exp(-D|\vec{q}|^{2}t)$ , where $D$ is the diffusion constant of the liquid (11). In a glass phase, however, where there is frozen (albeit random) structure, $S(\vec{q},t)$ in the thermodynamic limit remains finite as $t\rightarrow\infty$ , decaying only because of the finite size of the simulated system. We choose $|\vec{q}|=2\pi/b_{0}$ , with $b_{0}$ the bond length of the crystal structure, so that $S(\vec{q},t)$ will decay as rapidly as possible in the liquid phase. In addition to averaging over $t^{\prime}$ , we average over many orientations of $\vec{q}$ to get better statistics. The results are shown in Fig. 3, where the normalized time-dependent structure factor $S(|\vec{q}|,t)/S(|\vec{q}|,0)$ is plotted versus time delay $t$ , for four different temperatures around $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 $k_{B}T$ =0.40 and 0.45 eV, but remains finite for the lower temperatures, $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, $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 $k_{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 $T_{m}$ , and at $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 $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 $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=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 $F_{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.