## Abstract

We present a nonlinear model that allows exploration of the relationship between energy relaxation, thermal conductivity and the excess of low-frequency vibrational modes (LFVMs) that are present in glasses. The model is a chain of the Fermi–Pasta–Ulam (FPU) type, with nonlinear second neighbour springs added at random. We show that the time for relaxation is increased as LFVMs are removed, while the thermal conductivity diminishes. These results are important in order to understand the role of the cooling speed and thermal conductivity during glass transition. Also, the model provides evidence for the fundamental importance of LFVMs in the FPU problem.

## 1. Introduction

Glasses are solids that do not have long range order (Langer 2007). They are formed by cooling a melt fast enough to avoid crystallization. This process is not well understood (Jackle 1986; Anderson 1995; Phillips 1996). For example, it is not clear why the glass transition temperature (*T*_{g}), defined as the temperature at which the relaxation time exceeds the experimental time scale, can be lowered or raised by adding impurities (Micoulaut & Naumis 1999; Louzguine-Luzgin *et al.* 2008). Recently, it has been experimentally observed in metallic glasses that *T*_{g} is related to the thermal conductivity κ_{T}, which, at the same time, is related to the chemical composition (Louzguine-Luzgin *et al*. 2008). In fact, the glass-forming ability, determined by the minimal speed of cooling required to form a glass, depends on the chemical composition of the melt. Also, almost all glasses present anomalies in the density of low-frequency vibrational modes (LFVMs) (Elliot 1990; Binder & Kob 2005), such as the Boson peak (Binder & Kob 2005) or the floppy mode peak (Kamitakahara *et al.* 1991) owing to the flexible/rigid character of the atomic network (Phillips 1979; Thorpe 1983; Naumis 2005). It is surprising that most of the theories concerning glass transition do not give any special importance to such observations (Debenedetti 1996). In previous papers, we showed that *T*_{g} depends on the number of LFVMs, which is a function of the chemical composition due to the rigidity of the lattice (Naumis 2000, 2006; Huerta & Naumis 2002). It remains to be explained why the minimal cooling speed and the thermal conductivity are also related to the glass transition and chemical composition. This minimal cooling speed depends on the time for thermal relaxation. We must point out that, although a supercooled liquid is not a glass, it still shows a well-defined excess of LFVMs for time scales below structural rearrangements (Binder & Kob 2005). To study relaxation, nonlinear Hamiltonians (Fermi *et al.* 1955) such as the Fermi–Pasta–Ulam (FPU) model of a one-dimensional chain with nonlinear interaction are required (Fermi *et al.* 1955; Toda 1988; Lepri *et al.* 2003; Campbell *et al.* 2004). Still, there are many unsolved questions concerning this simple model. However, it seems that relaxation is dominated mainly by LFVMs (Ford 1961; Arnold & Moore 2001; Reigada *et al.* 2001; Ponno 2005), owing to a cascade energy transfer mechanism akin to turbulence. The purpose of this article is to answer some of these questions by looking at a simplified glass model, in which the number of LFVMs can be changed at will, mimicking what happens in a glass when the chemical composition is modified. The model is based on the application of rigidity theory (RT) to glasses, in which each covalent bond is considered to be a mechanical constraint (Naumis 2000). The present article also gives clues about the important role of LFVMs in the FPU model.

## 2. A nonlinear model of a glass

Our glass model is a one-dimensional lattice made of equal masses (*m*) joined by nonlinear springs, with second neighbour springs (SNSs) added at random as shown in figure 1. The corresponding Hamiltonian is
2.1
where *u*_{n} is the displacement of the mass *m* at site *n*, and *N* is the number of sites. *k* is the harmonic spring constant for first neighbours and *k*_{2} for second neighbours. and are the strengths of the nonlinear springs. In what follows *m*=1. We will use dimensionless variables for temperature, specific heat and heat flux. *Θ*_{n+2,n} is a random variable that takes values 0 or 1, with probabilities *c* and 1−*c*, respectively. Eventually, *Θ*_{n+2,n} can also be made a periodic function, as shown in figure 1*b*.

Let us briefly explain the characteristics of the model. When the anharmonic term is small, the actual number of LFVMs depends on the structure of the pure harmonic Hamiltonian, since, using the rotating phase approximation for small , one can show that the effect of nonlinearity is a shift of the modes (Cerón *et al.* 2005). Following the idea of RT, each bond can be treated as a mechanical constraint (Phillips 1979; Thorpe 1983). Comparing the number of constraints with the degrees of freedom, one can obtain the fraction of modes with zero frequency (known as floppy modes), denoted by *f*. In the case of a linear chain, removing bonds does not lead to a ‘floppy lattice’; it just cuts the lattice. So the only possibility is to add constraints, which is precisely the role of the springs that connect SNSs.

In the harmonic case (), the number of LFVMs changes as *c* goes from 0 to 1.0. There are many ways to prove this. For acoustic modes, ω(*q*)=*v*_{c}*q**a*, where *q* is the wavevector, *a* the lattice constant (set to 1 in the present article), ω(*q*) its corresponding frequency and *v*_{c} the speed of sound for a given *c*. Using that the density of states (DOS or ρ(ω)) in one dimension is given as ρ(ω)=(1/π)d*q*/dω(*q*), ρ(ω) turns out to be a constant, given by ρ(ω)=1/π*v*_{c}. figure 2 presents a plot of ρ(ω), and the corresponding limiting value 1/*v*_{c}π, obtained by diagonalizing the Hamiltonian and averaging over realizations of disorder, i.e. different configurations of SNSs. It can be seen in the inset that 1/*v*_{c} is a linear function that interpolates between ρ_{c=1}(ω) and ρ_{c=0}(ω).

## 3. Thermal conductivity

It is known that glass transition and chemical composition are related to the thermal conductivity (Louzguine-Luzgin *et al.* 2008). We can shed some light on the problem by looking at the effects of removing the LFVM. First, we treat the case of a pure harmonic Hamiltonian. The Kubo–Greenwood formula allows us to calculate the thermal conductivity (Elliot *et al.* 1974),
3.1
where is the phonon distribution function, Ω the system volume, *T* the temperature, *A*^{μ}≡1/2(*r*_{i}−*r*_{j})_{μ}Φ the transversal or longitudinal vibrational modes, Φ the dynamical matrix, μ and ν the Cartesian coordinates and *G*(ω)=(*m*ω^{2}−Φ)^{−1} the Green function. Such a function can be calculated by matrix inversion (Economou 1983). The harmonic model works only in the low-temperature region, because anharmonic effects are not present. These effects are the cause of the lattice thermal conductivity drop in dielectric materials after Umklapp processes activate.

figure 3 presents the thermal conductivity for chains with 100 atoms and different concentrations of SNSs using equation (3.1), in units where *k*_{B}=1 and . Note that the zero concentration case corresponds to the first neighbouring interaction. The second neighbouring interaction changes κ(*T*) with *c*, as can be seen in the inset of figure 3. Basically, κ(*T*) is reduced as *c* grows, until it reaches a minimum around *c*=0.5. Then κ(*T*) increases. This result can be explained as follows. κ(*T*) can be written as (Lepri *et al*. 2003; Binder & Kob 2005),
3.2
where *l*(ω) is the phonon mean free path, *C*(ω,*T*) the specific heat and *v*_{g}(ω) the group velocity. Let us analyse the limiting value at high temperatures, and, for the moment, we will consider a periodic arrangement of SNSs, as in figure 1*b*. In such a case, *C*(ω,*T*) =*k*_{B}ρ(ω)=(*k*_{B}/π)d*q*/dω(*q*). The group velocity is *v*_{g}(ω)=dω(*q*)/d*q*. Thus, in equation (3.2), the linear increase in the speed of sound as is compensated by a decrease in *C*(ω,*T*). As a result, equation (3.2) can be written as,
3.3
where ω_{S}(*c*) is the bandwidth of the phonon spectrum. For periodic chains, *l*(ω)=*N**a* and ω_{S}(*c*)≈(2*k*+*c**k*_{2}/2). This explains why the conductivity is higher for *c*=1 than for *c*=0, since ω_{S}(1) > ω_{S}(0) (figure 3). In the case of disordered chains, the arguments are similar except for two differences: *l*(ω) depends on the frequency and it is not possible to define a dispersion relation ω(*q*). However, the integral is dominated by extended low-frequency modes since *l*(ω)≈*N**a*. For LFVMs, the phonon spectrum is continuous, and it is possible to use an effective dispersion relation (Lepri *et al.* 2003). High-frequency modes are mainly localized and do not contribute much to the integral. Localized and acoustic modes are separated at the Ishii limit (Lepri *et al.* 2003) (denoted by ω_{I}(*c*)). For example, a detailed inspection of figure 2 for the concentration *c*=1/2 reveals that the Ishii limit is close to ω_{I}(1/2)≈1.4. Below this frequency, ρ(ω) is smooth owing to the extended nature of the excitations, while ρ(ω) is spiky above ω_{I}(1/2), a fact indicating localization. Therefore, κ(*T*) per unit length can be approximated by
3.4

One can see in figure 3 how κ(*T*) decreases as *c* goes from 0 to 1/2, as deduced from the fact that ω_{S}(0)>ω_{I}(1/2) (figure 2). Then it grows again after reaching a minimum. Such behaviour can be understood as a consequence of disorder, which is maximal at *c*=1/2, leading to strong localization. In real glasses and quasi-periodic systems at very low temperatures, the mean free path is also strongly reduced owing to Rayleigh scattering and two-level systems (Binder & Kob 2005). To study κ(*T*) for nonlinear cases, we solved the dynamical equations by using a fourth-order Runge–Kutta algorithm, with fixed boundary conditions (*u*_{0}=*u*_{N+1}=0). Both ends of a chain were in contact with heat baths. The equations obtained from equation (2.1) are (Lepri *et al.* 2003)
3.5
where δ_{nl} is a Kronecker delta. ξ_{+} and ξ_{−} are stochastic variables that follow a Langevin dynamics to simulate heat baths at temperature *T*_{+} and *T*_{−}, respectively. λ is the coupling constant between the bath and the chains (Lepri *et al.* 2003). The local heat flux is given by (Lepri *et al.* 2003)
3.6
and the total flux is . This leads to the thermal conductivity κ(*T*)=−*J*/∇*T*. Also, the local temperature is (Lepri *et al.* 2003)
3.7
where μ_{t} is the average over time. This allows us to build a temperature profile for the temperature gradient imposed on a chain. figure 4 presents such profiles for different concentrations of SNSs, using chains with and . Each profile has been averaged over 10 realizations of disorder. The profile shows that the chosen parameters are inside the Fourier law behaviour (Aoki & Kusnezov 2001). The inset presents the resulting κ(*T*) at *T*=1.0 as a function of *c*. The overall behaviour is similar to the linear case. These results allow us to understand in a qualitative way how chemical doping leads to changes in the thermal conductivity by affecting the LFVM.

## 4. Relaxation

The energy relaxation was studied by solving the motion equations with a fourth-order Runge–Kutta algorithm. First, the chains were thermalized using a Langevin dynamics, then the thermal bath was retired and a damping term was added at both ends of the chains. The equations of motion are
4.1
where Γ_{n}=γ[δ_{n,1}+δ_{n,N}] is the dissipation near the ends of the chain and γ is the damping. We used fixed-end boundary conditions. figure 5 presents a typical plot of the energy relaxation using a nonlinear Hamiltonian for different concentrations of SNSs. Each of the plots for the energy relaxation was made for *N*=100, starting from thermalized baths at *T*=0.5. An average of 50 realizations of disorder was made in all cases.

Clearly, the time required for relaxation increases as *c* goes from 0 to 0.5. A similar plot is obtained from *c*=0.5 to 1.0. The long time behaviour can be understood in terms of the depletion of the LFVM. According to Reigada *et al.* (2001), relaxation requires a transference of energy from high-frequency modes into LFVMs, which dissipate energy (Ponno 2005). The reduction of the LFVM as *c* increases means that not so many modes are available to dissipate energy and then energy relaxation is slower. Such a result is important in glasses because glass transition depends on a fast cooling that avoids thermal equilibrium. The chemical composition affects the number of LFVMs, and thus the speed required to achieve thermal equilibrium.

## 5. Conclusions

We presented a nonlinear model that shows in a clear way how low-frequency modes have a great impact on the thermal properties of a glass, although most of the theories concerning glass transition do not consider such important effects.

## Acknowledgements

We would like to thank C. Moukarzel for pointing out to us the rigidity properties of the second-neighbour interaction model. We thank DGAPA-UNAM projects IN-117806, IN-111906 and CONACyT 48783-F, 50368 and FICSAC-Universidad Iberoamericana for financial support.

## Footnotes

↵† On sabbatical leave from: Departamento de Física-Química, Instituto de Física, Universidad Nacional Autónoma de México (UNAM), Apartado Postal 20-364, 01000 México, Distrito Federal, Mexico.

One contribution of 14 to a Theme Issue ‘Topics on non-equilibrium statistical mechanics and nonlinear physics’.

- © 2009 The Royal Society