text:
# Vector vibrations and the Ioffe-Regel crossover in disordered lattices

## I Introduction

Vibrational excitations in disordered structures is an active field of investigation by different methods. The first approaches were developed many decades ago for the case of substitutional alloys (see e.g. Ref. Lifshitz\_64 ; Taylor\_67 ; Flicker\_73 ; Elliott\_74  and references therein). Disorder was introduced there by the random positioning of two or more atomic species (of different masses) onto lattice sites. Successful analytical methods have also been adapted for such problems Elliott\_74 ; Economou\_83:book ; Gonis\_92:book . The rigidity-percolation problem has been solved to some extent numerically and within the mean-field approximations Feng\_85 ; Garboczi\_85:VDOS . Atomic vibrations in structures with fractal geometries were investigated in Refs. Alexander\_82 ; Bunde\_92 . Some phenomenological Galperin\_89:review ; Klinger\_88:review ; Buchenau\_92 ; Heuer\_98  and more realistic Allen\_93 ; Feldman\_93 ; Fabian\_99  models of atomic vibrations in glasses have been developed in order to describe the thermodynamical and phonon-transport peculiarities of such materials. Various computer experiments have also been performed on disordered structures (see e.g. Taraskin\_99:PRB  and references therein).

Recently, the methods of the mode-coupling theory have been applied to investigate short-time vibrational dynamics, especially in the low-energy regime Gotze\_00 . The low-energy behaviour of the vibrational density of states and dynamical structure factor have been investigated in disordered lattice models within the coherent potential approximation (CPA) Schirmacher\_98 ; Martin\_Mayor\_00 ; Taraskin\_01:PRL , and in positionally disordered models treated perturbatively in the low inverse-particle-density limit Mezard\_99 ; Martin\_Mayor\_01:JCP ; Grigera\_01 ; Grigera\_01:BP . The spectral properties of the Euclidian random matrices associated with topologically disordered systems Mezard\_99  have been investigated within a supersymmetric statistic field theory in Ref. Offer\_00 . The spectral properties of the Laplacian defined on a random graph Biroli\_99  and ”small-world” lattices Monasson\_99  are related to topologically disordered solids and have also been investigated.

Previously, vibrational displacements in disordered structures have been treated as scalars, which is a convenient simplification of the vector nature of atomic vibrations (although vector vibrations were considered in the rigidity-percolation problem Feng\_85 ). However, a consistent treatment of vector vibrations in disordered lattices has not been performed up to now. The main purpose of this paper is to fill this gap and to give a mean-field description of vector vibrational models, namely: (i) to derive the CPA equations for vector vibrations in disordered lattices with force-constant disorder; (ii) to solve them analytically in limiting cases and numerically through the whole energy range; (iii) to apply the solution to investigate the scattering properties of external plane-wave vibrational excitations and, in particular, the Ioffe-Regel crossover from a weak- to strong-scattering regime; (iv) to demonstrate the importance of sum-rule correlations in the dynamical matrix for scattering properties of vibrational external plane waves. One of the main results of the paper is the phase diagram for vector vibrations in the representative f.c.c. lattice showing the ranges of weak- and strong-scattering regimes for different degrees of disorder. Another result concerns a comparative analysis of the vibrational problem with the electronic one, and an explanation of the essential differences in phase diagrams for these problems. All the major derivations and conclusions are supported by precise numerical solutions for the same model, which demonstrate the very good performance of the mean-field CPA approach, both for the density of states and the spectral density in a plane-wave basis (dynamical structure factor) in the whole energy range, except for the high-frequency tail region containing localized states.

The rest of the paper is arranged in the following manner. The model is defined in Sec. II. The properties of the effective mean field are analysed in Sec. III. The spectral densities and the scattering properties of the external plane waves are discussed in Sec. IV. The conclusions are given in Sec. V.

## II Model

### II.1 Hamiltonian

Atomic vibrations in the harmonic approximation can be treated in the classical limit Elliott\_74 . Thus the problem is reduced to the Hamiltonian formalism, in which the dynamical operator plays the role of the Hamiltonian, and the squared frequency, $`\omega ^2`$, is an energy, $`\epsilon =\omega ^2`$. The relevant Hamiltonian (dynamical operator) for vector vibrations of $`N`$ interacting atoms of unit masses ($`m_i=1`$) can be written in the following form:

$$\widehat{𝐇}=\underset{i\alpha ,j\alpha ^{}}{}H_{i\alpha ,j\alpha ^{}}|i,\alpha j,\alpha ^{}|,$$
(1)
where we have used the site orthonormal basis of vectors $`|i,\alpha `$, the only non-zero component of which is related to the displacement of atom $`i`$ ($`i=1,\mathrm{},N`$) along the Cartesian direction $`\alpha `$ ($`\alpha =1,\mathrm{},D`$). The (dynamical) matrix elements $`H_{i\alpha ,j\alpha ^{}}`$, being the second partial derivatives of the interatomic potential, obey the following sum rules Maradudin\_71 :

$$H_{i\alpha ,i\alpha ^{}}=\underset{ji}{}H_{i\alpha ,j\alpha ^{}}.$$
(2)
These sum-rule correlations between the elements of the $`D\times D`$ diagonal and off-diagonal blocks distinguish the vibrational problems from the electron problems described by similar Hamiltonians Economou\_83:book ; Gonis\_92:book  (see Sec. IV.5). The matrix elements of the Hamiltonian (1), generally speaking, are functions of all equilibrium atomic positional vectors, $`𝐫_i^{(0)}`$.

Our aim is to solve the eigenproblem for the Hamiltonian (1). This can be easily done for a crystal Maradudin\_71  but is not so simple for a disordered structure. The analysis starts with a choice of the vibrational model for a disordered system. This could be: (i) a crystalline lattice of atoms of random masses (mass disorder) and/or connected by random springs (force-constant disorder) Schirmacher\_98 ; Martin\_Mayor\_00 ; Taraskin\_01:PRL ; Kantelhardt\_01 ; (ii) positionally disordered structures created, for example, by random atomic displacements around crystalline sites Martin\_Mayor\_01:JCP ; Grigera\_01 ; Grigera\_01:BP , or by other ways, e.g. the Bethe lattice or by bond switching Barkema\_00 ; (iii) atomic structures, created e.g. by molecular dynamics (classical or ab-initio), in which intrinsic positional disorder defines the force-constant disorder (see e.g. Refs. Sarnthein\_95 ; Cappelletti\_95 ; Nakhmanson\_00 ; Vollmayr\_96 ; Taraskin\_97:PRB2 ). The simplest models are from the first class. Analytical approaches can be easily developed for them due to the existence of a well-defined reference crystalline system (i.e. the same lattice without disorder) Elliott\_74 ; Ehrenreich\_76 ; Economou\_83:book ; Gonis\_92:book . Surprisingly at first sight (see below), lattice models with force-constant disorder mimic very well the main features of topologically disordered glasses and even liquids, e.g. the boson peak, localized band tails, main peaks in the density of states and zero-energy singularity in the instantaneous spectrum Schirmacher\_98 ; Martin\_Mayor\_00 ; Taraskin\_01:PRL ; Simdyankin\_00 ; Taraskin\_02:PRB\_sing . The known drawback of lattice models is related to the $`k^4`$-frequency dependence of the width of the dynamical structure factor Sheng\_94 ; Schirmacher\_98 ; Martin\_Mayor\_00 ; Grigera\_01  as compared to the parabolic behaviour found experimentally (see e.g. Foret\_96 ). This $`k^4`$-dependence leads to a parabolic frequency scaling of the dynamic structure factor, $`S(k,\omega )\omega ^2`$ for $`\omega 0`$, which is in contrast to the frequency-independent $`S(k,\omega )`$ found for topologically disordered models Gotze\_00 ; Ruocco\_00 .

### II.2 Distribution of force constants

In this paper, we consider as a vibrational model the crystalline f.c.c. lattice of atoms, in which the nearest neighbours are connected by springs that are unstretched in equilibrium. The matrix elements of the Hamiltonian (1) in this case have the following form:

$$H_{i\alpha ,j\alpha ^{}}=\kappa _{ij}(\widehat{𝐫}_{ij}^{(0)})_\alpha (\widehat{𝐫}_{ij}^{(0)})_\alpha ^{},\text{if}ij,$$
(3)
where $`\widehat{𝐫}_{ij}^{(0)}=(𝐫_j^{(0)}𝐫_i^{(0)})/|𝐫_j^{(0)}𝐫_i^{(0)}|`$ is the unit connection vector between atoms $`i`$ and $`j`$ in equilibrium. The spring constants, $`\kappa _{ij}`$, are supposed to be independent of atomic positions and are random values taken from a certain probability distribution, $`\rho (\kappa )`$. Different functional forms of $`\rho (\kappa )`$ have been used for various problems, e.g. a combination of two $`\delta `$-functions Webman\_81 ; Odagaki\_81 ; Schirmacher\_99 ; Martin\_Mayor\_00 , (also used in the rigidity-percolation problem Feng\_85 ) and the normal (gaussian) Schirmacher\_98 , box Taraskin\_01:PRL  and hyperbolic Kantelhardt\_01  distributions for glassy problems. Below we use the box distribution:

$$\rho (\kappa )=\frac{1}{2\mathrm{\Delta }}\left[\theta (\kappa \kappa _0+\mathrm{\Delta })\theta (\kappa \kappa _0\mathrm{\Delta })\right],$$
(4)
where the Heavyside function $`\theta (x<0)=0`$ and $`\theta (x0)=1`$, $`\kappa _0`$ is the spring constant in the reference crystal and $`2\mathrm{\Delta }`$ is the full width of the distribution. The particular form of the peak-shaped distribution is not important for all results presented below from a qualitative point of view, as we have checked by using the normal distribution instead of the box one. The box distribution also has the advantage of not needing a lower cut-off to ensure positive values of force constants, as in the case for the normal and hyperbolic distributions.

If the distribution given by Eq. (4) is wide enough, i.e. $`\mathrm{\Delta }>\kappa _0`$ ($`\kappa _0=1`$ below), then negative spring constants appear in the system which could result in vibrational instability of the model in equilibrium, i.e. the occurrence of negative eigenvalues. Such a situation is of direct relevance Taraskin\_02:PRB\_sing  to an instantaneous mode analysis in liquids Wu\_92 ; Wan\_94 ; Keyes\_97:review .

The vibrational spectrum of a disordered system is defined by the properties of the Hamiltonian (dynamical) matrix. In the case of a crystal, the Hamiltonian matrix is well structured (e.g. cyclic band diagonal). For topologically disordered glasses with conserved local order, the Hamiltonian matrix is structured similarly to that for the crystalline counterpart (with the same local order) but now, due to the presence of topological disorder: (i) the positions of elements are distributed about their positions in the crystalline matrix (locational disorder) and (ii) the values of the elements are also distributed because of different interatomic distances. The first effect is completely ignored in the disordered-lattice models and the total influence of disorder is attributed entirely to the second effect. The question is, how good is this approximation?

One criterion relates to similarities in the probability distributions for the values of the matrix elements for both types of models. We have compared the distributions of the Hamiltonian matrix elements for two vibrational models: (i) the f.c.c. disordered-lattice model discussed above and described by the matrix elements $`H_{i\alpha ,j\alpha ^{}}^{\mathrm{lat}}`$ obeying Eq. (3) and (ii) a topologically disordered model of a single-component glass with predominantly icosahedral order (IC-glass) obtained by classical molecular dynamics (see Refs. Simdyankin\_00 ; Simdyankin\_01  for more detail). In the latter model, the Hamiltonian matrix elements, $`H_{i\alpha ,j\alpha ^{}}^{\mathrm{glass}}`$, are defined by the pairwise interatomic potential Dzugutov\_92:glass , $`V(|𝐫_{ij}^{(0)}|)`$, its spatial derivatives, $`V^{}(|𝐫_{ij}^{(0)}|)`$ and $`V^{\prime \prime }(|𝐫_{ij}^{(0)}|)`$, and the relative atomic positions,

$$H_{i\alpha ,j\alpha ^{}}^{\mathrm{glass}}=(\widehat{𝐫}_{ij}^{(0)})_\alpha (\widehat{𝐫}_{ij}^{(0)})_\alpha ^{}\left[V^{\prime \prime }(|𝐫_{ij}^{(0)}|)\frac{V^{}(|𝐫_{ij}^{(0)}|)}{|𝐫_{ij}^{(0)}|}\right]\delta _{\alpha ,\alpha ^{}}\frac{V^{}(|𝐫_{ij}^{(0)}|)}{|𝐫_{ij}^{(0)}|}.$$
(5)
The Hamiltonian matrix contains $`D\times D`$ diagonal and off-diagonal blocks, the diagonal ($`H_{i\alpha ,i\alpha }`$ and $`H_{i\alpha ,j\alpha }`$) and off-diagonal ($`H_{i\alpha ,i\alpha ^{}}`$ and $`H_{i\alpha ,j\alpha ^{}}`$) elements being distinctively distributed. For the diagonal elements in the off-diagonal blocks, it is convenient to compare the traces taken with opposite signs. They represent the force constants, $`\kappa _{ij}^{\mathrm{lat}}=_\alpha H_{i\alpha ,j\alpha }^{\mathrm{lat}}`$ and $`\kappa _{ij}^{\mathrm{glass}}=_\alpha H_{i\alpha ,j\alpha }^{\mathrm{glass}}=(V^{\prime \prime }(|𝐫_{ij}^{(0)}|)+2V^{}(|𝐫_{ij}^{(0)}|)/|𝐫_{ij}^{(0)}|)`$, for interactions beween atoms $`i`$ and $`j`$ in lattice and glassy models, respectively. Similarly the traces of diagonal blocks correspond to the total force constants, $`\kappa _{ii}^{\text{lat}}`$ and $`\kappa _{ii}^{\text{glass}}`$, for atom $`i`$ in lattice and glassy models, respectively, where $`\kappa _{ii}^{\text{lat(glass)}}=_j\kappa _{ij}^{\text{lat(glass)}}`$.

The four distributions for both models presented in Fig. 1(a)-(d) show qualitatively the same features. The distributions of the diagonal elements in the diagonal blocks, $`\rho (\kappa _{ii}^{\mathrm{glass}})`$ and $`\rho (\kappa _{ii}^{\mathrm{lat}})`$, are both peak-shaped with remarkably comparable peak positions and peak widths (see Fig. 1(a)). The distributions of the off-diagonal elements in the diagonal blocks are both centred around zero and have comparable widths (see Fig. 1(b)). The distributions of the elements for the off-diagonal blocks are much more model dependent. The solid line in Fig. 1(c) represents the distribution $`\rho (\kappa _{ij}^{\mathrm{glass}})`$ for all interacting atomic pairs. The distribution function has a complicated form, the features of which can be understood in terms of the functional forms of the interatomic potential and the atomic pair correlation function. For example, if nearest-neighbour interactions only are taken into account (as in the disordered-lattice model), then the distribution function $`\rho (H_{ij}^{\mathrm{glass}})`$ has a single peak (see the inset in Fig. 1(c)), which corresponds to the first peak in the pair correlation function Dzugutov\_92:glass . The simplest model distribution which can mimic the function $`\rho (\kappa _{ij}^{\mathrm{glass}})`$ (at least its scale) is the box distribution (the dashed line in Fig. 1(c)). The distribution of the off-diagonal elements in the off-diagonal blocks (see Fig. 1(d)) has a maximum at zero and is symmetric about it. The disordered-lattice models are intrinsically anisotropic, so that the distribution function, $`\rho (H_{i\alpha ,j\alpha ^{}}^{\mathrm{lat}})`$, can look quite different to $`\rho (H_{i\alpha ,j\alpha ^{}}^{\mathrm{glass}})`$. For example, for the ordered f.c.c. lattice ($`\mathrm{\Delta }=0`$), it consists of two $`\delta `$-functions symmetrically located around zero. With increasing disorder, these $`\delta `$-functions are transformed into two box distributions which merge at $`\mathrm{\Delta }=1`$, forming the single box distribution shown in Fig. 1(d) by the dashed line. For larger $`\mathrm{\Delta }>1`$, these two boxes overlap around zero, thus resulting in a box-shaped peak in this region, i.e. the distribution function $`\rho (H_{i\alpha ,j\alpha ^{}}^{\mathrm{lat}})`$ becomes closer in shape to the peak-shaped function $`\rho (H_{i\alpha ,j\alpha ^{}}^{\mathrm{glass}})`$.

The above comparison demonstrates that the very simple disordered-lattice model with a box distribution of force constants is able to reproduce correctly the ranges of the distributions of the Hamiltonian matrix elements for a topologically disordered glass with minimal short-range order (i.e. the IC-glass). This leads to the plausible conjecture that the physical origins of different phenomena, such as the boson peak, Ioffe-Regel crossover and secondary peak in the dynamic structure factor (analysed below), localization-delocalization transition, etc., are common for topologically disordered glasses and disordered lattices.

The other important point to emerge concerns the possibility of estimating the degree of disorder in disordered-lattice models at which they are able to mimic glassy models. For small disorder, all the distributions for the disordered-lattice model presented in Fig. 1 are quite narrow and are not comparable in width with the distributions for the glassy model. They become similar only for sufficiently large disorder ($`\mathrm{\Delta }1`$).

Below we use the simplest (box) distribution of force constants which allows us to analyse the problem analytically and which mimics some of the features of topologically disordered glasses by disordered lattices. A natural extension of the lattice models is in the use of more realistic force-constant distributions for crystalline counterparts of topologically disordered glasses and a comparison of the resulting properties of such models with the corresponding glasses. A possible way to mimic disorder in glasses is to introduce thermal (positional) disorder in the crystalline counterpart. We have done this for the IC-glass and its crystalline counterpart, the Frank-Kasper $`\sigma `$-phase Simdyankin\_02:PRB . A $`\sigma `$-phase lattice has been heated by means of molecular dynamics and the dynamical matrix for the instantaneous cofiguration has been calculated. The elements of such a matrix do not exhibit locational disorder and only their values are distributed. Using temperature as a control parameter it is possible to obtain the distributions of matrix elements very similar to those ones for the topologically disordered IC-glass (see Fig. 6 in Ref. Simdyankin\_02:PRB ). The vibrational density of states (VDOS) for instantaneous configurations of the $`\sigma `$-phase, as a consequence of such a similarity, resembles the VDOS for the IC-glass Simdyankin\_02:PRB ; Simdyankin\_00 . This supports the use of lattice models for mimicking topologically disordered structures. We should mention, however, that the elements of the dynamical matrix (belonging to the Euclidian ensemble Mezard\_99 ) for positionally disordered crystalline counterparts are correlated with each other within a certain row. These correlations are absent for the simple box distributions used below.

### II.3 Coherent potential approximation

The solution of the eigenproblem for the Hamiltonian (1) is known for a crystal, when $`\mathrm{\Delta }=0`$ in Eq. (4). The eigenmodes then are plane waves $`|𝐤,\beta `$, characterized by the wavevector $`𝐤`$, the branch number $`\beta `$ and the eigenenergies $`\epsilon _{𝐤\beta }`$. The VDOS, $`g^{\mathrm{cryst}}(\epsilon )=(DN)^1_{𝐤,\beta }\delta (\epsilon \epsilon _{𝐤\beta })`$, forms a band (see the dashed line in Fig. 2) containing van Hove singularities, and behaves as $`(\epsilon )^{(D/2)1}`$ around the band edges Maradudin\_71 . Disorder ($`\mathrm{\Delta }>0`$) brings new features to the spectrum (see Fig. 2): (i) the van Hove singularities are washed out; (ii) a high-frequency tail containing localized states occurs; (iii) extra states appear in the low-energy regime.

The exact solution of the vibrational problem with random spring constants is not known, but a number of approximate methods, both for electrons and vibrations, have been developed Elliott\_74 ; Economou\_83:book ; Efetov\_97:book ; Mirlin\_00:review . We are interested in the global features of the spectrum, except the tail region where the states are localized. An adequate approach in this case is a mean-field treatment within the coherent potential approximation (CPA) Elliott\_74 ; Ehrenreich\_76 ; Economou\_83:book , which we use below. The CPA is known to be very successful in a description of the electronic and vibrational properties of substitutional alloys and is expected to work very well in the whole energy range, except the far tails containing strongly localized states, where field-theoretical approaches are much better Mirlin\_00:review .

The Hamiltonian (1) can be rewritten in the bond representation,

$$\widehat{𝐇}=\frac{1}{2}\underset{i,j}{}\kappa _{ij}|ijij|,$$
(6)
with $`|ij=_\alpha (\widehat{𝐫}_{ij}^{(0)})_\alpha (|i,\alpha |j,\alpha )`$, which is a convenient form for a subsequent mean-field treatment. The VDOS of disordered systems, $`g(\epsilon )`$,

$$g(\epsilon )=\frac{1}{ND}\mathrm{Tr}\widehat{𝐀}(\epsilon )=\frac{1}{ND}\underset{𝐤\beta }{}\underset{d}{}|d|𝐤,\beta |^2\delta (\epsilon \epsilon _d),$$
(7)
can be defined via the spectral-density operator Ehrenreich\_76 , $`\widehat{𝐀}(\epsilon )=\delta (\epsilon \widehat{𝐇})`$, where the angular brackets, $`\mathrm{}`$, denote averaging over random spring constants. The coefficients, $`|d|𝐤,\beta |^2`$, in Eq. (7) have a simple physical meaning Ehrenreich\_76 ; they define the weight of the disordered eigenstate $`|d`$ in the crystalline state $`|𝐤,\beta `$, which is a plane wave of a certain polarization. The importance of the spectral-density operator in the plane-wave basis is related to the fact that the matrix elements $`\widehat{𝐀}_{𝐤\beta }(\epsilon )=𝐤\beta |\widehat{𝐀}|𝐤\beta `$ are proportional to the dynamical structure factor Martin\_Mayor\_00  (see Sec. IV.4) which can be measured experimentally Elliott\_90:book .

Both the VDOS, $`g(\epsilon )`$, and the spectral density, $`\widehat{𝐀}_{𝐤\beta }(\epsilon )`$, can be found for vector vibrations within the mean-field approximation (see Appendix A):

$$g(\epsilon )=\frac{1}{\pi }\text{Im}\left[\frac{1}{z(\epsilon )}G^{\mathrm{cryst}}\left(\frac{\epsilon }{z(\epsilon )}\right)\right],$$
(8)
and

$$A_{𝐤\beta }(\epsilon )=\frac{1}{\pi }\frac{z^{\prime \prime }(\epsilon )\epsilon _{𝐤\beta }}{\left[\epsilon z^{}(\epsilon )\epsilon _{𝐤\beta }\right]^2+\left[z^{\prime \prime }(\epsilon )\epsilon _{𝐤\beta }\right]^2},$$
(9)
via the crystalline Green’s function, $`G^{\mathrm{cryst}}`$, and the complex effective mean field, $`z(\epsilon )=z^{}(\epsilon )+iz^{\prime \prime }(\epsilon )=\stackrel{~}{\kappa }(\epsilon )/\kappa _0`$ (the dimensionless effective force constant), which can be found from the solution of the self-consistent Eq. (29).

## III Effective field

As seen from Eqs. (8)-(9), the spectral properties of a disordered lattice, within the CPA, are characterized by the energy dependence of the effective dimensionless force constant (effective field), $`z(\epsilon )`$. In this section, we discuss the properties of the effective field obtained by the solution of the self-consistent equation (29) (see Appendix A) for different values of disorder and in different energy ranges. The most important regime from the viewpoint of comparison with experimental data for the dynamical structure factor is the low-energy limit ($`\epsilon 0`$), analyzed below.

For the particular choice of the box-like force-constant probability distribution given by Eq. (4), the integration in Eq. (29) can be done analytically, resulting in the following expression:

$$1+\frac{1}{2\mathrm{\Delta }\alpha (\epsilon )}\text{ln}\left[\frac{1+\mathrm{\Delta }z(\epsilon )\alpha ^1(\epsilon )}{1\mathrm{\Delta }z(\epsilon )\alpha ^1(\epsilon )}\right]=0,$$
(10)
with

$$\alpha (\epsilon )\kappa _0ij|\stackrel{~}{\widehat{𝐆}}|ij.$$
(11)
Eq. (10) should be solved with respect to the complex effective field $`z(\epsilon )`$. Without disorder ($`\mathrm{\Delta }=0`$), the solution is trivial: $`z^{}=1`$ and $`z^{\prime \prime }=0`$. For finite disorder, the solution can be found numerically in the general case and analytically in some limiting cases.

### III.1 General case

The numerical solution of Eq. (10) for the real and imaginary parts of the effective field for different values of disorder in the whole energy range is presented in Fig. 3. As seen from Fig. 3(a), the real part of the effective spring constant varies around its crystalline value $`z^{}=1`$, being less than unity in the lower part of the energy band and greater than unity in the upper part of the band. This is the expected behaviour, because the value of $`z^{}`$ describes the level-repelling effect for the bare crystalline states when disorder is introduced in the system Taraskin\_01:PRL  (see Sec. IV.2). The deviation of $`z^{}`$ from unity increases with increasing disorder, which reflects the more pronounced degree of level repelling in more highly disordered systems.

The imaginary part of the effective field is non-zero and negative only in the band region. It basically reproduces the shape of the VDOS (Eq. (8)) and describes the widths of the spectral-density peak (Eq. (9)) which is proportional to the inverse life-time of plane waves in disordered structures (see Sec. IV.3). The magnitude of $`z^{\prime \prime }`$ increases with increasing disorder, thus indicating the shortening of the plane-wave life-time in strongly disordered lattices. The imaginary part of the effective field approaches zero in the low-energy limit, $`z^{\prime \prime }(\epsilon 0)0`$, only if the disorder is low enough, for $`\mathrm{\Delta }\mathrm{\Delta }_{}`$, when the structure in equilibrium is mechanically stable within the CPA (i.e. there are no negative eigenvalues).

### III.2 Low-energy limit

Using the energy dependence of the effective field in the low-energy regime (see Appendix B), the expression (8) for the VDOS can be rewritten as

$$g(\epsilon )=\frac{1}{z^{}(\epsilon )}g^{\mathrm{cryst}}\left(\frac{\epsilon }{z^{}(\epsilon )}\right)\frac{1}{\left(z^{}(0)\right)^{D/2}}g^{\mathrm{cryst}}(\epsilon )\frac{\chi _{\mathrm{Deb}}}{\left(z^{}(0)\right)^{D/2}}\epsilon ^{(D/2)1}.$$
(12)
As follows from Eq. (12), the low-energy VDOS in systems with lattice disorder functionally behaves the same as in the reference crystal, i.e. $`g^{\mathrm{cryst}}(\epsilon )\epsilon ^{(D/2)1}`$ (see Fig. 4), and differs from the crystalline VDOS just by the factor $`(z^{}(0))^{D/2}>1`$. In other words, extra states, in addition to the crystalline ones, appear in the low-energy regime due to disorder-induced level-repelling effects Ehrenreich\_76 ; Taraskin\_99:JPCM ; Taraskin\_99:PRB ; Taraskin\_01:PRL ; Finkemeier\_01 ; Simdyankin\_02:PRB . The relative density of these extra states is:

$$\frac{g(\epsilon )g^{\mathrm{cryst}}(\epsilon )}{g^{\mathrm{cryst}}(\epsilon )}\left[z^{}(0)\right]^{3/2}1.$$
(13)
With increasing energy, the disordered VDOS increases in energy faster than $`\epsilon ^{(D/2)1}`$ and this results in the boson peak, the origin of which has been discussed in Ref. Taraskin\_01:PRL .

In the low-disorder limit (see Appendix B), $`\mathrm{\Delta }0`$, the extra density of states (see Eq. (13)) in the low-energy range scales quadratically with disorder,

$$\frac{g(\epsilon )g^{\mathrm{cryst}}(\epsilon )}{g^{\mathrm{cryst}}(\epsilon )}\frac{D}{Z}\mathrm{\Delta }^2,\text{for}\mathrm{\Delta }1,$$
(14)
which has been confirmed numerically (see Fig. 5).

### III.3 Critical disorder

All the derivations presented above are valid only for sufficiently small disorder, $`\mathrm{\Delta }<\mathrm{\Delta }_{}`$. This restriction comes from Eq. (41) (see Appendix B), which has solutions only for $`\mathrm{\Delta }\mathrm{\Delta }_{}`$ (see Fig. 6). The critical value of disorder can be found numerically, giving $`\mathrm{\Delta }_{}1.296`$. The corresponding value of the real part of the effective field at criticality, $`z_{}^{}=z^{}(\mathrm{\Delta }_{})`$, can be found analytically using the fact that $`(df/dz^{})_{\mathrm{\Delta }=\mathrm{\Delta }_{}}=0`$ ( see Eq. (41)), so that

$$z_{}^{}=1\sqrt{1\frac{1\mathrm{\Delta }_{}^2}{\gamma }}.$$
(15)
Eq. (15) gives $`z_{}^{}0.43`$ for the f.c.c. lattice (with $`\gamma (Z/2D)1=1`$ for $`Z=12`$ and $`D=3`$).

When the disorder approaches the critical value where the VDOS behaves as a renormalized Debye function (see Eq. (12)) and the imaginary part of the effective force constant behaves as $`z^{\prime \prime }\epsilon ^{D/2}`$ (see Eq. (42)), the low-energy range shrinks to zero (see Fig. 4). Above the critical disorder, both the imaginary part of the effective field and the effective VDOS are finite at zero energy, $`z^{\prime \prime }(0)0`$ and $`g(\epsilon =0)0`$, and negative eigenvalues appear within the CPA. This means that the system becomes vibrationally unstable in equilibrium, even within the mean-field description. Such a situation corresponds to the spectrum of the dynamical matrix for instantaneous metastable configurations in liquid and glassy states. It turns out that the energy spectrum of the instantaneous dynamical matrix in the unstable regime ($`\mathrm{\Delta }>\mathrm{\Delta }_{}`$) exhibits a peculiar (singular) universal behaviour around zero energy, which we have also found in numerical experiments on liquid and vitreous SiO<sub>2</sub> and the IC-liquid and IC-glass Taraskin\_02:PRB\_sing .

## IV Spectral density

### IV.1 Shape of the spectral density

The spectral density in the plane-wave representation, $`A_{𝐤\beta }(\epsilon )`$, provides information about the distribution of contributions of the different disordered eigenstates in a particular plane wave from the branch $`\beta `$ and with wavevector $`𝐤`$. In a crystal, the spectral density is obviously a $`\delta `$-function, $`A_{𝐤\beta }^{\mathrm{cryst}}(\epsilon )=\delta (\epsilon \epsilon _{𝐤\beta })`$. Disorder introduced in the system transforms $`A_{𝐤\beta }^{\mathrm{cryst}}(\epsilon )`$ into the disordered spectral density which has a shape given by Eq. (9). The shape of $`A_{𝐤\beta }(\epsilon )`$ depends on the energy behaviour of the effective field $`z(\epsilon )`$. If we assume that the effective field is energy independent, then the spectral density has the shape of a Lorentzian peak, the width of which is proportional to the imaginary part, $`z^{\prime \prime }`$, of the effective field. Of course, the effective field does depend on energy (see Fig. 3), but if the imaginary part of the effective field is small enough, which is true in the whole energy range for sufficiently small disorder, and at least in the low-frequency regime ($`\epsilon _{𝐤\beta }\epsilon _{𝐤\beta }^{\mathrm{IR}}`$, see below) at any disorder below the critical one ($`\mathrm{\Delta }<\mathrm{\Delta }_{}`$), then the disordered spectral density still has the shape of a narrow peak (see Fig. 7). The position, $`\epsilon _{𝐤\beta }^{\mathrm{max}}`$, and full width at half-maximum, $`\mathrm{\Gamma }_{𝐤\beta }`$, of such a spectral-density peak can be estimated from the following equations:

$$\epsilon _{𝐤\beta }^{\mathrm{max}}z^{}(\epsilon _{𝐤\beta }^{\mathrm{max}})\epsilon _{𝐤\beta }\text{and}\mathrm{\Gamma }_{𝐤\beta }2|z^{\prime \prime }(\epsilon _{𝐤\beta }^{\mathrm{max}})|\epsilon _{𝐤\beta }.$$
(16)
For a weak energy dependence of the real part of the effective field (see Fig. 3(a)), we can approximate Eq. (16) giving:

$$\epsilon _{𝐤\beta }^{\mathrm{max}}z^{}(\epsilon _{𝐤\beta })\epsilon _{𝐤\beta }.$$
(17)
In the low-energy regime ($`\epsilon _{𝐤\beta }0`$), the imaginary part of the effective field is $`z^{\prime \prime }\epsilon ^{D/2}`$ according to Eq. (42), so that

$$\mathrm{\Gamma }_{𝐤\beta }\epsilon _{𝐤\beta }^{(D/2)+1}.$$
(18)
This equation corresponds to the well-known Rayleigh law for the peak width, $`\mathrm{\Gamma }_\omega (\omega _{𝐤\beta })`$, in the $`\omega `$-representation Sheng\_94 ; Schirmacher\_98 ; Martin\_Mayor\_00 :

$$\mathrm{\Gamma }_\omega (\omega _{𝐤\beta })=\frac{1}{2\omega }\mathrm{\Gamma }_{𝐤\beta }(\epsilon _{𝐤\beta }^2)\omega _{𝐤\beta }^{D+1},$$
(19)
or $`\mathrm{\Gamma }_\omega (\omega _{𝐤\beta })\omega _{𝐤\beta }^4`$ in the three-dimensional case (see Fig. 8).

Eq. (17) defines the dispersion law in disordered lattices, $`\omega _{𝐤\beta }^{\mathrm{max}}=\sqrt{\epsilon _{𝐤\beta }^{\mathrm{max}}}c_\beta ^{\mathrm{dis}}k`$ (for $`\epsilon _{𝐤\beta }0`$), which is basically the same as in the reference crystal except for the renormalized sound velocity, $`c_\beta ^{\mathrm{dis}}\sqrt{z^{}(0)}c_\beta ^{\mathrm{cryst}}`$, with $`c_\beta ^{\mathrm{cryst}}`$ being the sound velocity for branch $`\beta `$ in the reference crystal. Bearing in mind that $`z^{}(0)<1`$ (Fig. 3(a)), it is evident that $`c_\beta ^{\mathrm{dis}}<c_\beta ^{\mathrm{cryst}}`$ (see the inset in Fig. 8). We have found a similar relation for the model of the topologically disordered IC-glass, where $`c^{\mathrm{dis}}/c^{\mathrm{cryst}}0.72`$ Simdyankin\_02:PRB . This value corresponds to $`\mathrm{\Delta }1.28`$ (see the inset in Fig. 8), for which the distribution of the diagonal elements in the diagonal blocks for the disordered-lattice model fits very well the similar distribution for the IC-glass (see Fig. 1(a)).

### IV.2 Level-repelling effects

As seen from Eq. (17), disorder results in a shift of the crystalline spectral density, characterized by the real part of the effective field, $`z^{}`$. The bare level moves downwards (upwards) if $`z^{}<1`$ ($`z^{}>1`$), i.e. when it is located in the lower (upper) part of the band (see Fig. 3(a)). The shift of the bare level, $`\epsilon _{𝐤\beta }`$, by disorder is a consequence of level-repelling effects Ehrenreich\_76 ; Taraskin\_99:JPCM ; Taraskin\_99:PRB ; Taraskin\_01:PRL . Indeed, disorder introduced in the system results in the appearance of non-zero interaction matrix elements between the bare states. These interactions lead to standard quantum repelling between levels Landau\_63 ; Maradudin\_71 . If the interaction matrix elements are approximately the same for all levels, then the bare crystalline level, $`\epsilon _{𝐤\beta }`$, from the lower part of the energy band is shifted downwards just because the number of the states above the bare level (and repelling it downwards) is larger than the number of states below it (and repelling it upwards). This explains qualitatively why the value of the real part of the effective field is $`z^{}<1`$ in the lower part of the band and vice versa, $`z^{}>1`$ in the upper half of the energy band (see Fig. 3(a)).

We have proved quantitatively (analytically and numerically) such a picture of level-repelling effects in the VDOS for the force-constant-disordered f.c.c. lattice Taraskin\_01:PRL  and in the topologically disordered IC-glass Simdyankin\_01 . A similar effect has also been found numerically in models of disordered Si Finkemeier\_01 .

An example of the spectral densities for different energies of external plane waves at relatively high disorder ($`\mathrm{\Delta }=1`$) is shown in Fig. 7. Their widths at sufficiently large $`\epsilon _{𝐤\beta }`$ become comparable to their positions and plane waves at such energies are very short-lived quasiparticles Taraskin\_00:PRB\_IR1 ; Taraskin\_00:PRB\_IR2 . The level-repelling effects are clearly seen in pronounced shifts (to lower energies) of the peak positions with respect to $`\epsilon _{𝐤\beta }`$.

We have also presented in Fig. 7 the results of precise numerical solutions for the spectral densities obtained by the the kernel polynomial method (KPM) Silver\_97 . The very good agreement between the CPA and KPM results supports the reliability of the mean-field treatment in obtaining the spectral densities.

### IV.3 Ioffe-Regel crossover and localization

Disorder also broadens a bare level, so that the $`\delta `$-functions making up the crystalline VDOS are broadened into peaks of finite width, related to the imaginary part of the effective field (see Eq. (16)). This occurs because a plane wave $`|𝐤,\beta `$ is not an eigenstate in the disordered lattice. So it decays with time and can be treated as a quasiparticle having a finite life-time. This life-time is inversely proportional to the width of the peak of the spectral density and this determines the physical meaning of $`z^{\prime \prime }(\epsilon )`$ Ehrenreich\_76 .

The scattering properties of external plane waves in disordered lattices strongly depend on their bare energy. If the energy is low enough, the spectral-density peak is narrow compared to its position, and the quasiparticles live for a long time, thus indicating a weak-scattering regime. For higher energies, the width of the spectral density becomes comparable to the peak position (see Fig. 8). This means that the quasiparticles there decay quickly and the bare plane waves are in the region of strong scattering. The crossover between weak- and strong-scattering regimes is called the Ioffe-Regel crossover Ioffe\_60 ; Taraskin\_00:PRB\_IR1 ; Taraskin\_00:PRB\_IR2 . The energy of the bare plane wave, $`\epsilon _{𝐤\beta }^{\mathrm{IR}}`$, at which such a crossover occurs is the Ioffe-Regel crossover energy.

The quantitative condition for the Ioffe-Regel crossover is Taraskin\_00:PRB\_IR2 ,

$$\frac{\omega _{𝐤\beta }^{\mathrm{max}}}{2\pi }\tau _{𝐤\beta }1,$$
(20)
with $`\omega _{𝐤\beta }^{\mathrm{max}}\sqrt{\epsilon _{𝐤\beta }^{\mathrm{max}}}`$ standing for the frequency and $`\tau _{𝐤\beta }`$ being the life-time of the quasiparticle. The quasiparticle life-time can be easily evaluated in the weak-scattering regime from the time dependence of the probability, $`|𝐤,\beta ;t|𝐤,\beta ;0|^2`$, to find the system at time $`t`$ in the same plane-wave state $`|𝐤,\beta `$ as it was at $`t=0`$ Ehrenreich\_76 , i.e. $`|𝐤,\beta ;t|𝐤,\beta ;0|^2\mathrm{cos}^2(\omega _{𝐤\beta }^{\mathrm{max}}t)\mathrm{exp}(t/\tau _{𝐤\beta })`$, where

$$\tau _{𝐤\beta }=\frac{1}{|z^{\prime \prime }(\epsilon _{𝐤\beta })|}\sqrt{\frac{z^{}(\epsilon _{𝐤\beta })}{\epsilon _{𝐤\beta }}}\epsilon _{𝐤\beta }^{(D+1)/2}.$$
(21)
The inverse life-time, $`\tau _{𝐤\beta }^1`$, of course, coincides with the full width, $`\mathrm{\Gamma }_\omega (\epsilon _{𝐤\beta })`$, of the spectral density in the $`\omega `$-representation Taraskin\_99:PRB , $`\tau _{𝐤\beta }^1=\mathrm{\Gamma }_\omega (\epsilon _{𝐤\beta })\mathrm{\Gamma }_{𝐤\beta }/(2\sqrt{\epsilon _{𝐤\beta }^{\mathrm{max}}})`$, and the Ioffe-Regel criterion (20) can therefore be rewritten as

$$\frac{\omega _{𝐤\beta }^{\mathrm{max}}(\epsilon _{𝐤\beta })}{2\pi }\mathrm{\Gamma }_\omega (\epsilon _{𝐤\beta }).$$
(22)
This equation should be solved with respect to $`\epsilon _{𝐤\beta }`$ (or $`\omega _{𝐤\beta }=\sqrt{\epsilon _{𝐤\beta }}`$) and an estimate, $`\epsilon _{𝐤\beta }^{\mathrm{IR}}`$ (or $`\omega _{𝐤\beta }^{\mathrm{IR}}`$), for the Ioffe-Regel crossover energy (or frequency) can be found. Such a solution can exist at sufficiently large $`\omega _{𝐤\beta }`$. Indeed, in the low-frequency limit, the width of the spectral density, $`\mathrm{\Gamma }_\omega \omega _{𝐤\beta }^{D+1}`$, is much smaller than the peak position, $`\omega _{𝐤\beta }^{\mathrm{max}}\omega _{𝐤\beta }`$, but $`\mathrm{\Gamma }_\omega `$ grows quickly with increasing $`\omega _{𝐤\beta }`$ and can easily reach the peak-position value, if the disorder is large enough (see Fig. 8).

We have solved Eq. (22) numerically using the exact energy dependence for the effective field for different values of disorder, $`\mathrm{\Delta }`$ (see Fig. 8). The results for the Ioffe-Regel crossover energy are presented in the phase diagram for vector vibrations in the f.c.c. lattice (see Fig. 9). There is a minimum value of disorder ($`\mathrm{\Delta }_{\mathrm{min}}0.7`$ for the disordered f.c.c. lattice using Eq. (22)), starting from which the region of the strong-scattering regime appears in the system. The strong-scattering regime becomes broader with increasing disorder, while the low-energy weak-scattering regime shrinks basically to zero around the critical disorder, $`\mathrm{\Delta }_{}`$. We should stress that the boundary between weak- and strong-scattering regimes should really be thought of as a relatively wide crossover strip, and Eq. (22) serves for an order-of-magnitude estimate of the location of this strip. Bearing this in mind, the high-frequency region at sufficiently large disorder should not properly be regarded as a region of weak-scattering because the width of the spectral-density peak there is of the same order of magnitude as the peak position. The true weak-scattering regime occurs only either in the low-energy regime at any disorder less than the critical value, $`\mathrm{\Delta }\mathrm{\Delta }_{}`$, or for all energies at small disorder ($`\mathrm{\Delta }\mathrm{\Delta }_{\mathrm{min}}`$). The existence in these regions of small values of the energy and disorder parameters forces the spectral-density width to be much smaller than the peak position. In order to illustrate how sensitive the Ioffe-Regel crossover energy is to the numerical coefficients in Eq. (22), we have calculated $`\epsilon _{𝐤\beta }^{\mathrm{IR}}`$ from Eq. (22) in which the spectral-density peak width, $`\mathrm{\Gamma }_\omega `$, has been replaced, say, by $`2\mathrm{\Gamma }_\omega `$, thus enhancing the strong scattering. The results for this Ioffe-Regel crossover boundary are shown by the dot-dashed curve in Fig. 9. It is clearly seen how the region of the strong-scattering regime increases, with just the low-energy and low-disorder regions surviving for the weak-scattering regime.

We have also shown in Fig. 9 the trajectory of the boson peak (circled solid line) calculated as in Ref. Taraskin\_01:PRL . Remarkably, it almost coincides with the trajectory of the Ioffe-Regel crossover for large enough disorder when the strong-scattering regime appears in the system. This feature has also been found in real glasses, e.g. v-SiO<sub>2</sub> Taraskin\_00:PRB\_IR2  and in the topologically disordered IC-glass model (to be published elsewhere).

The solid line in the phase diagram shows the upper band edge of the VDOS calculated within the CPA. This appears to be a very good estimate for the trajectory of the threshold for the localization-delocalization transition, which, in fact, we have calculated Taraskin\_02:PhilMag  more exactly by means of a multi-fractal analysis Janssen\_98:review  with a relatively high precision, $`5\%`$ to be situated at slightly lower energies than the CPA band edge. The CPA fails to reproduce the high-energy tail in the VDOS (see Fig. 2) and more accurate methods should be used Economou\_83:book ; Gonis\_92:book . It should also be mentioned that, for $`\mathrm{\Delta }>1`$, a low-energy tail of localized states extending into the negative eigenvalue range appears in the true spectrum Taraskin\_02:PRB\_sing . The CPA is not able to reproduce this effect for $`1\mathrm{\Delta }\mathrm{\Delta }_{}`$ but still, for the majority of states with energies far enough away from both localization thresholds, the CPA results are very reliable.

### IV.4 Spectral density for large momentum transfer

In the previous subsections, we have analysed the situation at small and intermediate energies of external plane waves. These energies correspond to a relatively small momentum transfer in scattering experiments (see e.g. Foret\_96 ; Benassi\_96 ). However, larger values of $`k`$ are also of interest both from experimental (see e.g. Wischnewski\_98 ; Otomo\_98 ; Arai\_99 ; Nakamura\_01 ) and theoretical (see e.g. Ribeiro\_98 ; Horbach\_98:JNCS ; Taraskin\_00:PRB\_IR1 ; Taraskin\_00:PRB\_IR2 ; Gotze\_00 ; Ruocco\_00 ; Simdyankin\_01 ) points of view. Attention has been mainly focused on the behaviour of the dynamic structure factor, $`S_{𝐤\beta }(\omega )`$, which can be roughly approximated in the high-temperature regime ($`\mathrm{}\omega _\mathrm{d}/T1`$) by the following expression Ruocco\_00 :

$$S_{𝐤\beta }(\omega )\frac{Tk^2}{2m\mathrm{}\omega ^2}\underset{\mathrm{d}}{}|d|𝐤\beta |^2\delta (\omega \omega _\mathrm{d}),$$
(23)
where the Debye-Waller factor has been ignored (cf. Taraskin\_97:PRB1 ). In Eq. (23), we have also neglected the mixture of polarizations for phonons $`|𝐤\beta `$ and assume that $`\beta `$ is related to the longitudinal branch. In this case, the dynamical structure factor is obviously (see Eq. (7)) proportional to the spectral density,

$$S_{𝐤\beta }(\omega )\frac{Tk^2}{2m\mathrm{}\omega ^2}A_{𝐤\beta }^{(\omega )}\frac{Tk^2}{m\mathrm{}\omega }A_{𝐤\beta }(\omega ^2),$$
(24)
where $`A_{𝐤\beta }^{(\omega )}`$ is the spectral density in the $`\omega `$-representation.

A peculiar feature has been found in the frequency dependence of the dynamical structure factor in topologically disordered glasses at large values of the momentum transfer. Namely, the appearance of a secondary peak below the Brillouin peak in the $`\omega `$-dependence of $`S_{𝐤\beta }(\omega )`$ has been observed Horbach\_98:JNCS ; Gotze\_00 ; Ruocco\_00 ; Simdyankin\_01  with increasing value of $`k`$. This has been related to the boson peak Horbach\_98:JNCS , or interpreted as a manifistation of a microscopic relaxation process typical of a topologically disordered glass Ruocco\_00 .

A natural question is whether or not the secondary peak exists for the disordered-lattice models and, if so, then what is the origin of this peak. In order to answer these questions, first, we plot in Fig. 10(a) the evolution of the spectral density $`A_{𝐤\beta }(\epsilon )`$ for the longitudinal branch with increasing value of the wavevector, taken for definitness to be along the $`[100]`$ direction. It is clearly seen from this figure that the spectral density has a pronounced peak (Brillouin peak) for all values of $`k`$ taken from the first Brillouin zone. However this peak is no longer Lorentzian for large $`k`$. It is asymmetric and has a pronounced shoulder in the low-energy range. A similar picture has been observed in simulations for topologically disordered systems. The Brillouin peak shifts to higher energies, becomes very broad and a smooth shoulder-like feature appears below this peak and may also contain broad peak-shaped features Sampoli\_97 ; Ribeiro\_98 ; Taraskin\_00:PRB\_IR1 ; Simdyankin\_01 , so that the spectral density starts to resemble the VDOS. In the low-energy regime, the spectral density goes to zero.

The dynamical structure factor differs from the spectral density and contains an additional decaying energy (frequency) function (see Eq. (24). Therefore, if we plot the effective dynamic structure factor, $`\stackrel{~}{S}_{𝐤\beta }(\omega )=(2m\mathrm{}/T)S_{𝐤\beta }(\omega )=(k^2/\omega ^2)A_{𝐤\beta }^{(\omega )}=(k^2/2\sqrt{\epsilon })A_{𝐤\beta }(\epsilon )`$, versus frequency (see Fig. 10(b)), the low-frequency features are enhanced due to the $`\epsilon ^{1/2}`$ contribution. The low-energy shoulder in the spectral density of the disordered-lattice models transforms, for large $`k`$ values close to the Brillouin-zone boundary, to the secondary peak in the dynamic structure factor (see the solid line in Fig. 10(b)). The intensity of the secondary peak depends on the degree of disorder, so that it becomes very pronounced around the critical disorder (see Fig. 11(b)), although it does not appear for the spectral density in the energy representation (cf Fig. 11(a)). Therefore, the presence of the secondary peak in $`S_{𝐤\beta }(\omega )`$ of disordered-lattice models for large values of wavevector and large degrees of disorder is a distinctive feature of such models. Consequently, the occurrence of a similar peak in topologically disordered structures can hardly be attributed solely to topological disorder (cf. Ref. Ruocco\_00 ).

The origin of the secondary peak in the dynamic structure factor of disordered-lattice models can be understood from the qualitative analysis of Eqs. (9) and (24). In the low-frequency limit, $`\omega 0`$ (or $`\epsilon 0`$), the energy dependence of the spectral density in Eq. (9) is dictated by the imaginary part of the effective field, $`A_{𝐤\beta }|z^{\prime \prime }(\epsilon )|\epsilon ^{3/2}=\omega ^3`$, so that $`S_{𝐤\beta }A_{𝐤\beta }/\omega \omega ^2`$. The low-energy limiting behaviour of $`z^{\prime \prime }(\epsilon )`$ deviates approximately around the energy corresponding to the shifted (due to the level-repelling effect) lowest van Hove singularity, i.e. the first minimum in the energy dependence of $`z^{}(\epsilon )`$ (cf. Figs. 3(a) and (b)). Above this energy, the energy dependence of the imaginary part of the effective field is appreciably weaker and basically reproduces the shape of the VDOS (see Fig. 3(b)). This results in a weaker dependence of the dynamic structure factor, which shows, in this intermediate frequency range, either the shoulder-like behaviour, if the wavevector magnitude is not large enough and the low-frequency tail of the Brillouin peak plays a significant role (see the long and short dashed lines in Fig. 10), or a secondary peak for large values of $`k`$ (see the solid line in Fig. 10), when the influence of the Brillouin peak is suppressed by the $`\epsilon ^{1/2}`$ function in Eq. (24). Therefore, the position of the secondary peak in $`S_{𝐤\beta }(\omega )`$ for the disordered lattice models approximately coincides with the shifted lowest van Hove singularity and thus with the position of the boson peak Taraskin\_01:PRL . However, we should stress that the presence of the pronounced secondary peak in the dynamic structure factor for large values of $`k`$ does not mean a major contribution of plane waves with such $`k`$ to the boson peak, which is proportional to the sum of all plane waves (see Ref. Taraskin\_01:PRL  for more detail).

The $`\omega ^2`$-dependence of the dynamic structure factor in the low-frequency regime found for the disordered-lattice models is consistent with a $`k^4`$-dependence of the spectral density width and is in contrast to the constant (non-zero), limit of $`S_{𝐤\beta }(\omega )`$ found for topologically disordered glasses Gotze\_00 ; Ruocco\_00 ; Simdyankin\_01 , and thus can be considered as a drawback of lattice models.

### IV.5 Sum-rule correlations and plane-wave scattering

We have seen in the previous subsection that, independent of the degree of disorder, low-energy external plane waves are weakly scattered in a force-constant disordered lattice. This is a consequence of the essentially different dependences on the plane-wave energy of the peak width ($`\mathrm{\Gamma }_{𝐤\beta }\epsilon _{𝐤\beta }^{5/2}`$) and peak position ($`\epsilon _{𝐤\beta }^{\mathrm{max}}\epsilon _{𝐤\beta }`$) of the spectral density (see Eqs. (17,18)), which gives rise to the occurrence of weak scattering ($`\mathrm{\Gamma }_{𝐤\beta }\epsilon _{𝐤\beta }^{\mathrm{max}}`$) for $`\epsilon _{𝐤\beta }0`$. The question is how generic is this feature.

In order to answer this question, let us consider a more general formulation of the problem using the following Hamiltonian in the site representation:

$$\widehat{𝐇}=\underset{i}{}(\epsilon _i+\gamma \underset{ji}{}\kappa _{ij})|ii|\underset{i,ji}{}\kappa _{ij}|ij|.$$
(25)
The diagonal elements of the Hamiltonian have two constituents: random on-site ”energies”, $`\epsilon _i`$, and the correlation contribution proportional to the sum of random off-diagonal elements, $`\kappa _{ij}`$, in the same row. The parameter $`\gamma `$ controls the sum-rule correlations between diagonal and off-diagonal elements. If $`\gamma =1`$ and $`\epsilon _i=0`$, the Hamiltonian (25) describes the scalar vibrational problem, while the case $`\gamma =0`$ corresponds to the well-known electron Anderson Hamiltonian with on-site and off-diagonal disorder Economou\_83:book . The spectral properties of a similar Hamiltonian, but defined on a positionally random set of sites with deterministic transfer integrals, have been studied in Ref. Mezard\_99 .

The scattering properties of external plane waves in the disordered structure described by the Hamiltonian (25) can be studied using spectral densities for which the mean-field approach works very well (as shown above for the vibrational problem). It is easy to show that only in the case of $`\gamma =1`$ and $`\epsilon _i=0`$ can the Hamiltonian (25) be reduced to the form of Eq. (6) which has a multiplicative effective mean field (the energy of the quasiparticles is the product of the effective field and the bare energy, $`\stackrel{~}{\epsilon }_{𝐤\beta }=z(\epsilon )\epsilon _{𝐤\beta }`$ \- see Eq. (32)). For $`\gamma 1`$, the situation is quite different, because now the effective field has an additive constituent as well, and this drastically changes the scattering properties of plane waves. Below, we demonstrate this, considering the simplest particular case for the Hamiltonian (25) characterized by $`\gamma =0`$, $`\kappa _{ij}=\kappa _0`$ and random $`\epsilon _i`$ taken from the box distribution of half-width, $`\mathrm{\Delta }^{\mathrm{el}}`$, centred around $`\epsilon _0=0`$, i.e. the standard Anderson Hamiltonian with diagonal disorder Kramer\_93 , for the simple cubic lattice.

The CPA spectral density, $`A_𝐤^{\mathrm{el}}(\epsilon )=(\pi )^1\mathrm{Im}(\epsilon \stackrel{~}{\epsilon }_𝐤)^1`$, for the Anderson Hamiltonian can be expressed via the quasiparticle energy, $`\stackrel{~}{\epsilon }_𝐤=\epsilon _𝐤+\stackrel{~}{\epsilon }(\epsilon )`$, which is now the sum of the crystalline energy, $`\epsilon _𝐤`$, and the effective field, $`\stackrel{~}{\epsilon }(\epsilon )=\stackrel{~}{\epsilon }^{}(\epsilon )+i\stackrel{~}{\epsilon }^{\prime \prime }(\epsilon )`$. The spectral densities $`A_𝐤^{\mathrm{el}}(\epsilon )`$ are peak shaped and in the low-energy limit, $`\delta \epsilon _𝐤=(\epsilon _𝐤\epsilon _{\mathrm{min}}^{\mathrm{cryst}})0`$ (the energy $`\delta \epsilon _𝐤`$ is referred to the lower crystalline band edge, $`\epsilon _{\mathrm{min}}^{\mathrm{cryst}}=\kappa _0Z`$, where $`Z`$ is the coordination number), the position and width of the peak can be estimated from the following relations:

$$\delta \epsilon _{\mathrm{max}}^{\mathrm{el}}\delta \epsilon _𝐤,$$
(26)
and

$$\mathrm{\Gamma }^{\mathrm{el}}(\epsilon _𝐤)\frac{1}{2}\pi \chi ^{\mathrm{el}}(\mathrm{\Delta }^{\mathrm{el}})^2\sqrt{\delta \epsilon _{\mathrm{max}}^{\mathrm{el}}}(\mathrm{\Delta }^{\mathrm{el}})^2\sqrt{\delta \epsilon _𝐤},$$
(27)
where the peak position, $`\delta \epsilon _{\mathrm{max}}^{\mathrm{el}}=(\epsilon _{\mathrm{max}}^{\mathrm{el}}\epsilon _{\mathrm{min}}^{\text{band}})`$, is measured from the CPA lower band boundary, $`\epsilon _{\mathrm{min}}^{\mathrm{band}}`$ (see Fig. 12). The last expression is a consequence of the fact that $`\text{Im}\left[G^{\mathrm{cryst}}(\epsilon \stackrel{~}{\epsilon })\right]\pi \chi ^{\mathrm{el}}\sqrt{\epsilon \epsilon _{\mathrm{min}}^{\mathrm{cryst}}}`$, if $`\left(\epsilon \epsilon _{\mathrm{min}}^{\mathrm{cryst}}\right)0`$ with $`\chi ^{\mathrm{el}}`$ being the constant coefficient similar to the coefficient $`\chi _{\mathrm{Deb}}`$ in Eq. (37).

In full analogy to vibrations (see Eq. (22)), the Ioffe-Regel crossover between weak- and strong-scattering regimes for electrons can be defined in terms of the following condition:

$$\frac{\delta \epsilon _{\mathrm{max}}^{\mathrm{el}}(\epsilon _𝐤)}{2\pi }\mathrm{\Gamma }^{\mathrm{el}}(\epsilon _𝐤).$$
(28)
In contrast to vibrations, the spectral-density peak width for electrons scales more slowly than the peak position as a function of the plane-wave energy (see Eqs. (26)-(27) and Fig. 12). This means that, in the long-wavelength limit, $`\delta \epsilon _𝐤0`$, the width of the spectral density is always larger than the peak position, and the strong-scattering regime occurs for plane waves having energies around the lower crystalline band edge (see Fig. 13). The situation is opposite for vibrations, for which $`\mathrm{\Gamma }_{𝐤\beta }\epsilon _{𝐤\beta }^{5/2}`$ (cf. Eqs. (18) and (27)). As a consequence, the shapes of the regions associated with the weak- and strong-scattering regimes for plane waves in the phase diagram are different for electrons and vibrations (cf. Figs. 9 and 13). Now, for electron plane waves with energies around the lower crystalline band edge, the scattering is always strong and the weak-scattering regime is realized only in the midband range for sufficiently small disorder. In contrast to vibrations, only a single parameter, the strength of disorder, $`\mathrm{\Delta }^{\mathrm{el}}`$, rather than two parameters, the degree of disorder $`\mathrm{\Delta }`$ and the energy, $`\epsilon _{𝐤\beta }`$, for vibrations, controls the appearance of the weak-scattering regime. The boundary between the weak- and strong-scattering regimes (see the dashed line in Fig. 13) is again not precise and should be thought of as an estimate for the crossover energy.

It should be mentioned that the electron phase diagram presented in Fig. 13, for the case of very small disorder ($`\mathrm{\Delta }^{\mathrm{el}}0`$), supports a standard concept for electrons Mott\_79:book , according to which the Ioffe-Regel crossover from the weak- to the strong-scattering regime corresponds to the transition from extended states to localized ones. Indeed, the Ioffe-Regel crossover boundary (the dashed line in Fig. 13) almost coincides with the localization-delocalization threshold (the circles in Fig. 13) Bulka\_87  and the CPA band edge (solid line) as $`\mathrm{\Delta }^{\mathrm{el}}0`$. At high disorder, the localization-delocalization transition, however, is not associated with the Ioffe-Regel crossover, because all plane waves for such disorder are strongly scattered.

In the more general case of diagonal and off-diagonal disorder coexisting in the Anderson Hamiltonian (but with $`\gamma =0`$), the most successful homomorphic cluster CPA Yonezawa\_79 ; Li\_88  results in an effective field which necessarily has an additive constituent and the above conclusions still hold. In the case of $`0<\gamma <1`$, the mean field can be constructed as well but it also has an additive constituent which disappears only for $`\gamma =1`$ (to be published elsewhere). Thus, we can conclude that the existence of the low-energy weak-scattering regime is characteristic only of the vibrational problem and it results from the exact sum-rule correlations in the corresponding Hamiltonian.

## V Conclusions

We have analytically investigated classical harmonic atomic vibrations in disordered lattices. The f.c.c. atomic lattice with force-constant disorder described by a box probability distribution has been chosen for analysis. The distributions of matrix elements in diagonal and off-diagonal blocks in the dynamical matrix are rather similar to those of topologically disordered models. The atomic vibrations were treated as vectors and analysed within the framework of the single-bond coherent potential approximation (CPA); the results are essentially the same as those found from a precise numerical study. The CPA self-consistency equation has been derived for vector models and solved analytically in the limiting cases of small disorder and low energies.

The CPA solution of the problem has allowed us to investigate the spectral properties of the model in terms of spectral densities. We have shown that the Rayleigh law for the spectral-density width is intrinsic for vector vibrations in disordered lattices.

We have also investigated the Ioffe-Regel crossover between weak and strong scattering for vector vibrations. The regime of weak scattering takes place at all energies for small degrees of disorder, but in the low-energy regime occurs only for large disorder. The existence of a disorder-independent weak-scattering regime for low-energy external plane waves is a consequence of the sum-rule correlations in the dynamical matrix, typical for the vibrational problem only. We have demonstrated that the Ioffe-Regel crossover occurs in the boson-peak region and is not related to the localization-delocalization transition for vector vibrations in the disordered lattice studied. In contrast, the Ioffe-Regel crossover for electrons (described by the Anderson model with on-site energy disorder in the simple cubic lattice) exists only at low values of disorder, where it is associated with the localization-delocalization transition. At higher values of disorder for the Anderson model, the transition is between strongly-scattered extended states and localized states.

We have found that a secondary peak appears in the dynamic structure factor below the Brillouin peak for sufficiently large values of $`k`$ for the lattice models as previously found in simulations of topologically disordered glasses. Hence this feature cannot be due solely to the presence of topological disorder, nor can it be due to relaxations.

## Appendix A

Here we give the derivation of the self-consistency CPA equation for vector vibrations. It is known Webman\_81 ; Odagaki\_81  that the self-consistency equation for the effective force constant, $`\stackrel{~}{\kappa }(\epsilon )`$, in the single-bond approximation can be written as follows:

$$\frac{\delta \kappa }{1\delta \kappa ij|\stackrel{~}{\widehat{𝐆}}|ij}|ijij|=0,$$
(29)
where $`\delta \kappa =\kappa _{ij}\stackrel{~}{\kappa }`$ and the effective Green’s function, $`\stackrel{~}{\widehat{𝐆}}=\left[\epsilon \stackrel{~}{\widehat{𝐇}}\right]^1`$, is defined for the effective Hamiltonian, $`\stackrel{~}{\widehat{𝐇}}=(\stackrel{~}{\kappa }/2)\underset{i,j}{}|ijij|`$. Our aim is to find the matrix element, $`ij|\stackrel{~}{\widehat{𝐆}}|ij`$, for vector vibrations. This can be rewritten as

$$ij|\stackrel{~}{\widehat{𝐆}}|ij=2\underset{\alpha ,\alpha ^{}}{}(\widehat{𝐫}_{ij}^{(0)})_\alpha (\widehat{𝐫}_{ij}^{(0)})_\alpha ^{}\left(\stackrel{~}{G}_{i\alpha ,i\alpha ^{}}\stackrel{~}{G}_{i\alpha ,j\alpha ^{}}\right),$$
(30)
and this can then be expressed in terms of diagonal elements of the Green’s function in the site basis for the reference crystal. In order to show this, let us first rewrite the matrix elements, $`\stackrel{~}{G}_{i\alpha ,j\alpha ^{}}`$ in the plane-wave basis, $`\{|𝐤,\beta \}`$,

$$\stackrel{~}{G}_{i\alpha ,j\alpha ^{}}=\frac{1}{N}\underset{𝐤\beta }{}(\widehat{𝐧}_\beta )_\alpha (\widehat{𝐧}_\beta )_\alpha ^{}\frac{\mathrm{exp}\left\{i𝐤(𝐫_i^{(0)}𝐫_j^{(0)})\right\}}{\epsilon \stackrel{~}{\epsilon }_{𝐤\beta }},$$
(31)
where $`\widehat{𝐧}_\beta `$ is a $`D`$-dimensional unit polarization vector for the branch $`\beta `$ and the effective energy, $`\stackrel{~}{\epsilon }_{𝐤\beta }`$, obeys the relation:

$$\stackrel{~}{\epsilon }_{𝐤\beta }=z(\epsilon )\epsilon _{𝐤\beta }=\stackrel{~}{\kappa }\underset{j(i)}{}\left(𝐫_{ij}^{(0)}\widehat{𝐧}_\beta \right)^2\left(1\mathrm{exp}\{i\mathrm{𝐤𝐫}_{ij}^{(0)}\}\right),$$
(32)
with $`z(\epsilon )=\stackrel{~}{\kappa }(\epsilon )/\kappa _0`$ being the dimensionless effective force constant and the index $`j(i)`$ running over all nearest neighbours to the site $`i`$, i.e. $`j=1,\mathrm{}Z`$, with $`Z`$ standing for the number of nearest neighbours (all the sites in a crystal are equivalent and $`\epsilon _{𝐤\beta }`$, in fact, does not depend on $`i`$).

Substituting expression (31) back into Eq. (30) we obtain:

$$ij|\stackrel{~}{\widehat{𝐆}}|ij=\frac{2D}{Z\stackrel{~}{\kappa }}\left(\epsilon \stackrel{~}{G}(\epsilon )1\right),$$
(33)
where

$$\stackrel{~}{G}(\epsilon )\stackrel{~}{G}_{ii}(\epsilon )\frac{1}{D}\underset{\alpha }{}\stackrel{~}{G}_{i\alpha ,i\alpha }(\epsilon )=\frac{1}{ND}\underset{𝐤\beta }{}\frac{1}{\epsilon \stackrel{~}{\epsilon }_{𝐤\beta }}$$
(34)
is the average on-diagonal (independent of $`i`$) element of the effective Green’s function in the site basis. The effective Green’s function, $`\stackrel{~}{G}(\epsilon )`$, can be expressed via the crystalline one, $`G^{\mathrm{cryst}}(\epsilon )`$, as

$$\stackrel{~}{G}(\epsilon )=\frac{1}{z}G^{\mathrm{cryst}}(\epsilon /z).$$
(35)
Finally, the expression for $`ij|\stackrel{~}{\widehat{𝐆}}|ij`$ reads:

$$ij|\stackrel{~}{\widehat{𝐆}}|ij=\frac{2D}{Zz\kappa ^{(0)}}\left[\frac{\epsilon }{z}G^{\mathrm{cryst}}\left(\frac{\epsilon }{z}\right)1\right],$$
(36)
and this solves the problem when substituted back into Eq. (29).

The dimensionless effective spring constant $`z(\epsilon )=\stackrel{~}{\kappa }/\kappa _0`$ fully determines the effective mean field. This is a complex quantity and should be found from the solution of the self-consistent equation (29). The self-consistent equation (29) for vector vibrations is very similar to that for the scalar models Webman\_81 ; Odagaki\_81 ; Schirmacher\_98 ; Martin\_Mayor\_00 . The only difference is the factor $`D`$ in the denominator of Eq. (29), with $`ij|\stackrel{~}{\widehat{𝐆}}|ij`$ given by Eq. (36). For scalar models, this is absent, even for $`D`$-dimensional lattices with $`D1`$.

The crystalline Green’s function, $`G^{\mathrm{cryst}}(x)`$, with complex argument $`x=\epsilon /z`$, can be easily found as an analytical continuation from the real axis of its imaginary part (VDOS) Ehrenreich\_76 . After that, the self-consistent equation can be solved numerically in the general case, and analytically in some limiting cases (e.g. $`\epsilon 0`$ and/or $`\mathrm{\Delta }0`$), and the effective field can be found.

If $`z(\epsilon )`$ is known, the task of evaluating the spectral density becomes straightforward. Indeed, bearing in mind that the energy, $`\stackrel{~}{\epsilon }_{𝐤\beta }`$, of the eigenstates in the effective crystal is connected with the energy, $`\epsilon _{𝐤\beta }`$, of the similar eigenstates in the reference crystal according to relation (32), it is easy to see that the VDOS and diagonal elements of the spectral-density operator in the crystalline plane-wave basis satisfy Eqs. (8)-(9).

## Appendix B

Here we show how the self-consistent Eq. (10) can be solved analytically in the low-energy limit, $`\epsilon 0`$, and how the effective field, $`z(\epsilon )`$, can be found in this case. In order to do this, first we find the expression for the function $`\alpha (\epsilon )=\alpha ^{}(\epsilon )+i\alpha ^{\prime \prime }(\epsilon )`$ given by Eq. (11) via $`z(\epsilon )`$ and $`\epsilon `$. This can be done, bearing in mind the following relations for the crystalline Green’s function, $`G^{\mathrm{cryst}}(\epsilon /z)`$, in the low-energy limit, when the Debye approximation for the crystalline VDOS,

$$g^{\mathrm{cryst}}\chi _{\mathrm{Deb}}\epsilon ^{(D2)/2},$$
(37)
is valid:

$$\text{Re}\left[G^{\mathrm{cryst}}\left(\frac{\epsilon }{z}\right)\right]\text{const}\text{and}\text{Im}\left[G^{\mathrm{cryst}}\left(\frac{\epsilon }{z}\right)\right]\pi \chi _{\mathrm{Deb}}\left(\frac{\epsilon }{z^{}}\right)^{\frac{D2}{2}}.$$
(38)
Substituting Eq. (38) into Eqs. (11) and (36), we obtain the leading terms in $`\epsilon `$ for $`\alpha `$:

$$\alpha ^{}(\epsilon )\frac{2D}{Zz^{}},$$
(39)
and

$$\alpha ^{\prime \prime }(\epsilon )\frac{2D}{Z(z^{})^2}\left[\frac{\pi \chi _{\mathrm{Deb}}}{z^{}}\left(\frac{\epsilon }{z^{}}\right)^{\frac{D}{2}}z^{\prime \prime }\right].$$
(40)
These expressions can then be used for an expansion of the logarithm in Eq. (10) in terms of $`\epsilon `$, which, after some algebra, results in the following final relations for $`z^{}(\epsilon )`$ and $`z^{\prime \prime }(\epsilon )`$:

$$f(z^{})\frac{2\mathrm{\Delta }}{z^{}(\gamma +1)}\text{ln}\left[\frac{1+\mathrm{\Delta }+\gamma z^{}}{1\mathrm{\Delta }+\gamma z^{}}\right]=0,$$
(41)
and

$$z^{\prime \prime }(\epsilon )\frac{2\pi D\chi _{\mathrm{Deb}}}{Zz^{}}\left(\frac{\epsilon }{z^{}}\right)^{\frac{D}{2}}\left\{\left[\mathrm{\Delta }^2+\left((\gamma +1)z^{}\right)^2\left(1+\gamma z^{}\right)^2\right]^1\left[(\gamma +1)(z^{})^2\right]^1\right\}^1,$$
(42)
where $`\gamma (Z/2D)1`$.

In fact, relation (41) is the equation to find the real part of the effective force constant, $`z^{}`$. The important point is that all the parameters entering this equation are energy independent, so that $`z^{}(\epsilon 0)z^{}(0)=\text{const}`$. Eq. (41) can be solved numerically in the case of arbitrary disorder and analytically in the small-disorder limit.

The imaginary part of the effective force constant, as follows from Eq. (42), scales with energy as $`z^{\prime \prime }(\epsilon )\epsilon ^{D/2}`$. Such a behaviour is basically a consequence of the Debye law for the VDOS in the low-energy regime (see Eq. (37)) and, in this sense, is general for disordered lattices. This dependence also determines the scaling behaviour of the spectral-density peak width with energy, the well-known Rayleigh law in this case.

Eqs. (41)-(42) have important consequences for the energy dependence of the VDOS in the low-energy limit. Indeed, in the limit $`\epsilon 0`$, the expression (8) for the VDOS can be rewritten in the form of Eq. (12).

In the low-disorder limit, $`\mathrm{\Delta }0`$, Eq. (41) can be solved with respect to $`z^{}`$ while Eq. (42) can be simplified, so that

$$z^{}(\epsilon )1\frac{2D}{3Z}\mathrm{\Delta }^2(1+O(\epsilon )),$$
(43)
$$z^{\prime \prime }(\epsilon )\frac{2\pi \chi _{\mathrm{Deb}}D}{3Z}\epsilon ^{D/2}\mathrm{\Delta }^2.$$
(44)
These relations are used in deriving Eq. (14).

————————————————————————