# Transitions between metastable states in silica clusters

## I Introduction

The theoretical investigation of the properties of disordered systems is a very difficult task, because the lack of simmetry usually prevents analytical approaches. Many authors have found it convenient to rely on the concept of potential energy landscape in the configuration space (see e.g. Frauenfelder et al. 1991, Berry 1993, Heuer and Silbey 1993, Mohanty et al. 1994, Heuer 1997, Angelani et al. 1998, Mousseau and Barkema 1998, Wales et al 1998).

From a purely qualitative point of view, we may imagine the multidimensional analogue of a surface rich in climbing points (first order saddles) and valleys (the various minima). A detailed analysis of the energy landscape topology requires numerical simulation.

In this work we analyze the potential energy landscape of SiO<sub>2</sub> clusters, namely \[SiO<sub>2</sub>\]<sub>20</sub>, \[SiO<sub>2</sub>\]<sub>30</sub>, and \[SiO<sub>2</sub>\]<sub>50</sub>, which are comparable in size to (or larger than) monoatomic and binary systems so far investigated. After a short description of the numerical methods used to move up and down the multidimensional hypersurface, we report our results regarding:

1. the structure of the clusters (resembling that of the bulk solid, although affected by surface effects);
2. the energy distribution of the stationary points of the potential energy function;
3. a physical interpretation of geometrical quantities uselful to characterize transitions between different minima;
4. the identification of tunneling systems.

## II Numerical methods

We simulate the interaction among Si and O atoms with the recently developed pair potential by Van Beest et al. (1990), modified with a very short range contribution (Guissani and Guillot 1996), necessary to avoid unphysical divergencies. It is given by

$$\mathrm{\Phi }_{ij}=\frac{q_iq_j}{r_{ij}}+a_{ij}\mathrm{exp}(b_{ij}r_{ij})\frac{c_{ij}}{r_{ij}^6}+4ϵ_{ij}\left[\left(\frac{\sigma _{ij}}{r_{ij}}\right)^{24}\left(\frac{\sigma _{ij}}{r_{ij}}\right)^6\right],$$
(1)   
where both i and j indexes run on all Si and O atoms, and the values of the parameters were determined by Van Beest et al. (1990) and Guissani and Guillot (1996). This potential has already been widely used (with slight modifications) in molecular dynamics studies of the liquid and glassy phases of silica (see e.g. Vollmair et al. 1996, Horbach et al. 1996, Taraskin and Elliott 1997). As our purpose is to simulate clusters, we used free boundary conditions. We adopted the procedure described in detail in Daldoss et al. (1998, 1999) to locate minimum-saddle-minimum triplets, that are the key ingredients to describe jump processes. Schematically the steps are:

* descent towards a minimum by the conjugate gradient method (Press et al. 1986), starting from a randomly chosen configuration;
* ascent towards the vicinity of a saddle, following the eigenvector corresponding to the lowest eigenvalue;
* once the potential energy along the path of the previous item starts decreasing, we take the corresponding configuration (hill-climbing point, Berry 1993) as the starting point for a Newton-Raphson search (Press et al. 1986) of the first order saddle point;
* a descent on both sides of the saddle following again the eigenvector of the minimum eigenvalue.

This technique provides approximate adjacent minima that are subsequently fed into the Newton algorithm for accurate location. With this numerical procedure we found many thousand minimum-saddle-minimum triplets (hereafter double well potentials, DWP).

## III Results

### A The Structure of the clusters

We obtain information on the structure of the cluster by the radial distribution function $`g(r)`$, defined in the same way as for extended systems: it shows the bond lengths, the short range order, and also allows for the calculation of the co-ordination number. The results are reported in Fig.1, which was obtained by averaging over the various configurations obtained (both minima and saddles) to have good statistics. We have plotted the partial $`g(r)`$ for the different bonds (Si-Si, Si-O, O-O) and cluster sizes. The resulting bond lengths obtained by the main peaks are in good agreement with the experimental data on bulk structures (see table I) and with simulations with periodic boundary conditions (Taraskin and Elliott 1997). Nevertheless we observe smaller peaks on the low-distance side of the nearest-neighbour peak in the Si-Si and O-O bonds. In our opinion this is due to surface effects: in fact it appears from the figure that these anomalies tend to disappear with increasing system size. Very satisfactory is also the result concerning the co-ordination number: it is 4 for Si-O bond, in agreement with the tetrahedral structure typical of SiO<sub>2</sub>. This is also confirmed by the analysis of Si-Si, also equal to 4, and O-O, equal to 6.

### B Energy distributions of stationary points

We now analyze the topological features of the potential energy hypersurface; at first we present the energy distribution of the stationary points (saddles and minima) that form double wells. In Fig. 2 we report the three energy distributions for saddles, lower minima and upper minima of the various DWP, respectively; only the result for 150 atoms is presented, the situation being quite analogous in the other two cases. We note that the distributions are very similar, with defined peaks superimposed onto a broad backgroud; the peaks are due to the fact that certain configurations are favoured for given values of the binding energy, and so they act like attraction basins during the descent towards the minima. The presence of these structures in the distribution is more evident with increasing system size: in \[SiO<sub>2</sub>\]<sub>20</sub> the curve looks smoother (Brangian 1998). It should be remarked that Fig. 2 does not refer to the total distribution of stationary points, because their number increases exponentially with the system size: we sampled partially the configuration space, and we are quite confident not to have introduced systematic errors in this sampling. The only thing we can note is that maybe our numerical investigation is in some way prevented from reaching very low lying (that is crystalline-like) configurations: since we are not interested in a careful thermodynamical analysis (Doye and Wales 1998), but only in transitions between stable disordered states, we think this fact does not constitute a serious drawback.

### C Topological features of double well potentials

Double well potentials can be characterized by many quantities: the first one is the asymmetry $`\mathrm{\Delta }`$, that is the energy difference between the two connected minima of a DWP. This parameter is essential for the identification of candidate two level systems (TLS), which require a value of $`\mathrm{\Delta }`$ of the order of less than $``$1 K (however, as we will see, this condition is not sufficient to identify a TLS). Our results indicate that the asymmetries are distributed exponentially, the most part being lower than 5000 K; the distribution is not sensitive to the system size. Equally important are the energy barriers $`V`$, i.e. the energy differences between a minimum and the corresponding saddle. Of course for every DWP there are two values of $`V`$, one for the relaxation process and one for the activation. The distribution of the relaxation barriers similar in shape to that of the $`\mathrm{\Delta }`$’s, but on average the values of $`V`$ are smaller, being significantly present only up to 1000 K. There seems to be little correlation between the asymmetry and the barrier (both for the activation and relaxation).

We have evaluated also the (mass weighted) euclidean distance between pairs of minima, defined as

$$dist(a,b)=\left[\underset{i}{}\frac{m_i}{\stackrel{~}{m}}|𝐫_{i,a}𝐫_{i,b}|^2\right]^{1/2},$$
(2)   
where $`a`$ and $`b`$ are the two minima, $`m_i`$ the single atom mass, $`𝐫_{i,\alpha }`$ its position in the $`\alpha `$ configuration, and $`\stackrel{~}{m}=_im_i/N`$. In Fig. 3 we report the relative distributions. The distributions do not extend very much beyond $``$10$`\sigma `$ and present a maximum at $`2\sigma `$.

Another very interesting parameter to consider is the participation number defined as

$$N_{\mathrm{part}}=\underset{i}{}\frac{d_i^2}{d_{\mathrm{max}}^2}.$$
(3)   
Here $`d_i`$ refers to the atom-atom distances in (2) and $`d_{max}`$ is the distance of the atom that moves most during the transition; the partecipation number gives an idea on the numbers of atoms involved in a transition. In Fig. 4(a) we report this quantity for the three cases studied; the same quantity normalized to the number of particles constituting the clusters is reported in Fig. 4(b). We see that $`N_{\mathrm{part}}`$ is a nearly scaling quantity with $`N`$, indicating that at least an appreciable part of the atoms that move in the transitions belong to the bulk.

We can make use of multidimensional transition state theory (for a review see Hanggi 1985) to estabilsh a link between the potential energy landscape and the (classical) relaxation dynamics of the system. Under appropriate conditions (Hanggi 1985) the classical, thermally activated transition probability between two metastable minima is given by

$`\mathrm{\Gamma }=\nu ^{}\mathrm{exp}\left({\displaystyle \frac{E_b}{\kappa _BT}}\right)`$ (4)   
$`\nu ^{}={\displaystyle \frac{\underset{i=1}{\overset{N}{}}\nu _i^M}{_{j=1}^{N1}\nu _j^S}},`$ (5)   
where $`\nu ^{}`$ is an effective frequency that takes into account the effect of all the degrees of freedom on the one which forces the transition. This frequency can be rewritten as

$$\nu ^{}=\frac{\nu _0}{R}$$
(6)   
with $`\nu _0`$ the lowest frequency of the dynamical matrix in the starting minimum, while $`R`$ is the product of all the $`N1`$ positive eigenvalues of the dynamical matrix at the saddle point, divided by the corresponding product at the minimum. In Fig. 5 we show the value of this entropic factor $`R`$ (which enhances or reduces the transition probability), versus the energy barrier value, both for the relaxation and the activation transitions. This plot evidences a rather marked correlation between $`R`$ and the barrier height, and in most cases $`R<`$1 or even $``$1.

### D Two-level systems (TLS)

We describe now the most difficult part of our work, that is the selection of two-level systems. Following the energy landscape paradigm, these are DWP that imply purely quantum mechanical relaxation processes. The calculation has been carried out by assuming the validity of the 1D Wentzel-Kramers-Brillouin (WKB) approximation (Froman and Froman 1965, Landau and Lifchitz 1967; Schiff 1968; a discussion on the validity of this approximation in the case of Ar clusters, as well as on the effects to be expected when it is released, can be found in Daldoss et al. 1999): the splitting of the ground state is given in terms of the action integral between the two minima. In principle, in order to apply the WKB procedure it would be necessary to find the least action path, i.e. the classical path that takes from one minimum to the other and minimizes the action integral; this involves rather heavy numerical calculations. As a starting point, we decided to use the path that takes from one minimum to the other, and that is defined at each point by the direction of the minimum eigenvector, as described in Section II. Work on the minimization of the action integral is in progress (Brangian et al. 1999). It should be stressed that as cumbersome as this procedure may look, it is probably the simplest way to get a quantitative estimate of the tunneling splitting. It should also be mentioned that, in agreement with previous works (Heuer and Silbey 1993, Heuer 1997, Daldoss et al. 1998), we find roughly one TLS for 1000 DWP’s, which implies a very extensive search strategy.

We have performed an a posteriori test on the reliability of the use of the 1D WKB scheme. The main hypothesis is that, along the chosen 1D path (the least-action or a very close one), the degrees of freedom other than the considered one are decoupled from it: in this case it is reasonable to expect that the Schroedinger equation may be (nearly) factorized into 3$`N`$ mutually independent equations, a condition for the applicability of the WKB scheme in many dimensions (Schiff 1968). If this factorization should actually take place, the eigenvalues of the dynamical matrix (which are proportional to the curvatures of the hypersurface) other than the lowest one should remain constant along the chosen path. The departure of these eigenvalues from constancy gives a measure (though a qualitative one) of the invalidity of the 1D scheme. In Fig. 6 we show, as an example, the results for two TLS belonging to clusters of 90 atoms. In one case (Fig. 6(b)) the approximation is satisfied in a very good way; on the contrary, in the other case (Fig. 6(a)) there is appreciable mixing of the low energy eigenvectors. These results are in qualitative agreement with those relative to Ar clusters (Daldoss et al, 1999); they imply that the actual structure of TLS in disordered systems is probably much more complex than expected, since many-dimensional effects seem to play important roles.

## IV Conclusions

We have reported preliminary results of an extensive investigation on the properties of the potential energy landscape in SiO<sub>2</sub> clusters of three different sizes (60,90, and 150 atoms); the aim of this research is to find connections between the properties of the landscape and the high- and low-temperature relaxation dynamics. By analyzing the structure of the clusters and the topological features of their energy landscape (and in particular of its minima and first-order saddle points), we have identified tunneling centres and studied the conditions of validity of the WKB approximation, which allows a quantitative estimate of the tunneling splitting.

ACKNOWLEDGMENTS

This work was supported in part by the Parallel Computing Initiative of the INFM.

REFERENCES

Angelani, L., Parisi, G., Ruocco, G., and Viliani, G., 1998, Phys. Rev. Lett., 81, 4648.

Berry, R.S., 1993, Chem. Rev., 93, 2379.

Brangian, C., 1998, Thesis (University of Trento).

Brangian, C., Pilla, O., and Viliani, G., 1999, to be published.

Daldoss, G., Pilla, O., and Viliani, G., 1998, Phil. Mag. B, 77, 689.

Daldoss, G., Pilla, O., Viliani, G., Brangian, C., and Ruocco, G., 1999, to appear on Phys. Rev. B.

Doye, J.P.K., and Wales D.J., 1998 Phys. Rev. Lett., 80, 1357.

Frauenfelder, H., Sligar, S.G., and Wolynes, P.G., 1991, Science, 254, 1594.

Froman, N., and Froman, O.O., 1965, JWKB Approximation (North-Holland: Amsterdam).

Guissani, Y., and Guillott, B., 1996, J. Chem. Phys., 104, 7633.

Hanggi, P., 1986, J. Stat. Phys., 42, 105.

Heuer, A., 1997, Phys. Rev. Lett., 78, 4051.

Heuer, A., and Silbey, R.J., 1993, Phys. Rev. Lett., 70, 3911.

Horbach, J., Kob, W., Binder, K., and Angell, C.A., 1996, Phys. Rev. E, 54, R5827.

Landau, L., and Lifchitz, E., 1967, Mécanique Quantique, chapter 7 (MIR: Moscow).

Mohanty, U., Oppenheim, I., and Taubes, C.H., 1994, Science, 266, 425.

Mousseau, N., and Barkema, G.T., 1998, Phys. Rev. E, 57, 2419.

Press, W.H., Flannery, B.P., Teukolsky, S.A., and Vetterling, W.T., 1986, Numerical Recipes (Cambridge University Press: Cambridge).

Schiff, L.J., 1968, Quantum mechanics (McGraw-Hill: New York).

Taraskin, S.N., and Elliott, S.R., 1997, Phys. Rev. B, 56, 8605.

Van Beest, B.W.H., Kramer, G.J., and Van Santen, R.A., 1990, Phys. Rev. Lett., 64, 1955.

Vollmair, K., Kob, W., and Binder, K., 1996 Phys. Rev. B 54, 15808.

Wales, D.J., Miller, M.A., and Walsh, T.R., 1998, Nature, 394, 758.

FIGURE CAPTIONS   

Fig. 1. Radial distribution function $`g(r)`$ in clusters of three different sizes. Besides the three main peaks, which are also found in bulk systems, we notice other peaks for the O-O and Si-Si bonds, probably due to the presence of surface co-ordination defects.

Fig. 2. Energy distribution for saddles, upper and lower minima in clusters of 150 atoms.

Fig. 3. Distribution of euclidean distances between connected minima, in units of $`\AA `$ and $`\sigma =1.6\AA `$, i.e. the Si-O bond length.

Fig. 4. Participation number: the $`y`$ axis units are chosen such that the integrals of the various curves is equal to 1. The bottom plot refers to the participation ratio normalized to the number of atoms.

Fig. 5. Entropic ratio $`R`$ as a function of the barrier height for $`N=`$60, 90, and 150.

Fig. 6. Variation of the 10 lowest eigenvalues of the dynamical matrix along the minimum eigenvalue path (see text) in two TLS of clusters of 90 atoms. The minima are arbitrarily assigned the $`\pm `$1 values of the coordinate along the path.