# The dynamical structure factor in disordered systems

## Abstract

We study the spectral width as a function of the external momentum for the dynamical structure factor of a disordered harmonic solid, considered as a toy model for supercooled liquids and glasses. Both in the context of single-link coherent potential approximation and of a single-defect approximation, two different regimes are clearly identified: if the density of states at zero energy is zero, the Rayleigh $`p^4`$ law is recovered for small momentum. On the contrary, if the disorder induces a non vanishing density of states at zero energy, a linear behaviour is obtained. The dynamical structure factor is numerically calculated in lattices as large as $`96^3`$, and satisfactorily agrees with the analytical computations.

PACS numbers: 61.43.Fs,63.50.+x

The spectrum of vibrational excitations of supercooled liquids and glasses is attracting a great deal of attention from the experimental side , the numerical simulations  and also from the analytical point of view . Vibrational excitations on the GHz range develop in the supercooled liquid, by the same temperatures where anomalous behaviour of the specific-heat is found. Some rather substance-independent features in this vibrational spectrum are found. For example, the vibrational density of states presents an excess respect to the usual Debye behaviour ($`g(\omega )\omega ^2`$) known as the Boson Peak. The dynamical structure-factor, $`S(p,\omega )`$, reveals well-defined sound-like peaks for wave-lengths not much greater than the inter-particle distance. Moreover, the speed of this high-frequency sound is close to the one of the low-frequency typical one. For fixed external momentum, $`p`$, the width of the spectral peak grows as $`p^2`$, which has recently been recovered in the mode coupling approximation (see ). However, in Ref. , a simpler model of a disordered three-dimensional harmonic solid was studied. There it was claimed that if some spring-constants are allowed to have a small negative value, (but not so much that negative-energy modes appear) a boson peak develops. In the context of the CPA approximation , the usual $`p^4`$ width of $`S(p,\omega )`$ arisen from Rayleigh-scattering is reported. It was, however, noticed that a rather different scaling appeared at the characteristic frequencies of the Boson-Peak. Furthermore, in Ref. , it was shown a $`p^2`$ broadening for a one-dimensional disordered harmonic solid model. In this paper, we want to investigate the problem, both by numerical and analytical means. Our main finding will be that the Rayleigh-scattering $`p^4`$ broadening holds, unless the solid becomes unstable. That is, if the dynamical-matrix, to be defined later, has an extensive number of negative eigenvalues, the system do present sound-like peaks in its $`S(p,\omega )`$, but with a width proportional to $`p`$, at least for very small values of $`p`$. Otherwise, the standard $`p^4`$ behaviour is to be expected.

To be more specific, our dynamical matrix in $`D`$ dimensions is given by

$`H`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{xy}{}}\varphi _x_{xy}\varphi _y,`$ (1)   
$`_{xy}`$ $`=`$ $`{\displaystyle \underset{\mu =1}{\overset{D}{}}}{\displaystyle \frac{1+\alpha _{y,\mu }}{2}}\left(\delta _{xy}\delta _{x,y+\widehat{\mu }}\right)`$ (2)   
$`+`$ $`{\displaystyle \frac{1+\alpha _{y\widehat{\mu },\mu }}{2}}\left(\delta _{xy}\delta _{x,y\widehat{\mu }}\right),`$ (3)   
where $`\widehat{\mu }`$ is the lattice unit-vector in the $`\mu `$ direction. Notice that we cannot separate the longitudinal and transverse modes since, for the sake of simplicity, our matrix has no internal indices as in Ref. . The $`\alpha _{y,\mu }`$ is the random part of the spring constant that joins the sites $`y`$ and $`y+\widehat{\mu }`$. In a finite lattice, periodic boundary conditions are applied. The dynamical matrix, verifies the constraint related with traslational symmetry: the vector of all equal components is eigenvector with zero eigenvalue. For a harmonic solid, the object that we will study is the dynamical structure function, $`S(p,E)`$, in the energy rather than in the frequency-domain:

$$S(p,E)=\frac{1}{\pi }\mathrm{Im}\underset{ϵ0}{lim}\overline{p|\frac{1}{E+iϵ}|p}.$$
(4)   
In the above expression, the bra-ket notation has been used and $`|p`$ stands for a normalized plane-wave of momentum $`p`$. As usual, the overline represents the average over the random variables $`\alpha `$. We choose to work in the energy, rather than in the frequency domain, since we mainly want to consider the case where a significant fraction of the spectrum is negative, and so the transformation $`E=\omega ^2`$ is no longer well defined.

We study the case where the $`\alpha _{y,\mu }`$ are uncorrelated, random variables whose probability distribution is

$$p(\alpha )=(1\rho )\delta (\alpha )+\rho f(\alpha ),0\rho 1.$$
(5)   
In the above equation, $`\rho `$ is the probability of finding one defect, while $`f`$ is a continuous probability function, that we take flat between $`\lambda `$ and $`0`$ ($`\lambda <0`$). Notice that the disorder of Ref.  is recovered, by taking $`\rho =1`$ and $`f`$ Gaussian. The rationale for choosing this kind of disorder is that a spring of an unusually large negative spring-constant induces a negative energy eigenvalue of the dynamical matrix $``$. In fact, if one keeps in Eq.(3) only the spring connecting sites $`1`$ and $`2`$, corresponding to the large negative defect, the eigenvectors are easily shown to be $`\varphi _x=(\delta _{x,1}+\delta _{x,2})/\sqrt{2}`$ and $`(\delta _{x,1}\delta _{x,2})/\sqrt{2}`$. However, assuming that the surrounding springs have their ordered value, one finds that the Hamiltonian (1) is negative for the displacement configuration $`(\delta _{x,1}\delta _{x,2})/\sqrt{2}`$ if the spring-constant verifies $`\alpha <(D+\frac{1}{2})`$. Therefore, one has a contribution of order $`\rho `$ to the density of states over the negative spectrum. Moreover, if we assume that some of the surrounding spring have very small positive values, a contribution of higher-order in $`\rho `$ to the density of states is generated, no matter how small is the negative value of the spring constant, $`\alpha +1`$. Therefore, if the probability of $`\alpha <1.0`$ is non-vanishing, the hybridization of the localized field-configuration described above, with the plane-waves eigenstates will non-trivially modify the eigenvectors of $``$. To study this, we shall perform a single-defect calculation of the resolvent.

In this way, we learn that the order $`\rho `$ threshold for the presence of negative eigenvalues is not $`\alpha =D1/2`$, as roughly shown in the introduction, but $`D`$. This threshold separate two well defined scaling limits for the width of the $`S(q,E)`$ at small $`q`$:

* If the density of states is null at zero energy the imaginary part of the self-energy $`\mathrm{\Sigma }`$ is proportional to $`p^{D+2}`$ (Rayleigh scattering, that yields a spectral width $`p^{D+1}`$ in the frequency domain).
* When an extensive number of negative eigenvalues is present, $`\mathrm{\Sigma }p^2`$ (or $`p`$ for the width in the frequency domain). Of course, the $`p^{D+2}`$ contribution will be still present but sub-dominant at low momentum. A crossover might be visible, depending on the strength of the disorder.

The dynamical matrix in the presence of an unique defect of amplitude $`\alpha `$ can be split in two terms:

$$_{xy}=_{xy}^0+_{xy}$$
(6)   
where the pure crystal matrix $`^0`$ has plane waves eigenvectors with eigenvalues

$$E_0(p)=\underset{\nu }{}(1\mathrm{cos}p_\nu ).$$
(7)   
The perturbation, hence, connects the two sites $`y^0`$ and $`y^0+\nu `$:

$`_{xy}`$ $``$ $`\alpha ^{y^0;\nu }|y^0;\mu y^0;\mu |,`$ (8)   
$`|y^0;\mu `$ $``$ $`{\displaystyle \frac{\delta _{x,y^0}\delta _{x,y^0+\nu }}{\sqrt{2}}}.`$ (9)   
The propagator can be written as

$$\frac{1}{z}=\frac{1}{z^0}+\frac{1}{z^0}𝒯\frac{1}{z^0},$$
(10)   
where the resummation of the harmonic series gives

$$𝒯|y^0;\nu \frac{\alpha }{1\alpha a(z)}y^0;\nu |,$$
(11)   
with

$$a(z)\frac{1}{D}\frac{d^Dq}{(2\pi )^D}\frac{E_0(q)}{zE_0(q)}.$$
(12)   
Notice that the correction term in Eq.(11), has an isolated singularity for the value of $`z`$ satisfying $`1\alpha a(z)=0`$. This value decrease monotonically from $`z=0`$, for $`\alpha =D`$, and will be called $`z_\alpha `$. This isolated singularity, correspond to the ground state of the dynamical matrix, the residue being $`\mathrm{\Psi }_F(x)\mathrm{\Psi }_F^{}(y)`$. Therefore, one obtain the wave-function for this eigenvalue

$$\mathrm{\Psi }_F(x)\frac{d^Dq}{(2\pi )^D}\frac{1e^{iq_\nu }}{z_\alpha E_0(q)}e^{iqx},$$
(13)   
that has a localization length of order $`|z_\alpha |^{1/2}`$. We thus see, that unless $`\alpha `$ was exceedingly close to the critical value $`D`$, the eigenvector is strongly localized around the defect. The single-defect approximation to the self-energy, Eq.(14), amounts to consider that each defect only contributes to its own localized eigenvector, and that no other defect is within its localization length. Let us turn back to the original problem, with an extensive number of defects, and repeat the above calculation, neglecting all terms which contain two different defects. The matrix $`𝒯`$ is now a sum of terms like the one in Eq.(11). If we now perform the average over the $`\alpha `$’s and apply Dyson resummation, the self-energy part of the propagator turns out to be

$$\mathrm{\Sigma }(z,p)=\rho E_0(p)𝑑\alpha f(\alpha )\frac{\alpha }{1\alpha a(z)},$$
(14)   
at first order in $`\rho `$ (the defect interactions will generate the order $`\rho ^2`$ and higher order corrections to the single defect result).

The width of the $`S(q,E)`$ is simply given by the value of the imaginary part of the self-energy at the peak, whose position can be obtained from the real part of the self-energy, $`E^{\mathrm{max}}(p)E_0(p)+\mathrm{Re}\mathrm{\Sigma }(E_0(p),p)`$. Notice that our self energy is proportional to $`E_0(p)`$ so that we are basically getting a finite renormalization of the speed of sound. In order to estimate the imaginary part of the self energy, it is useful to realize that for small positive values of the energy, one has

$$a(E+iϵ)\frac{1}{D}i\frac{E^{D/2}}{2^{D/2}\pi ^{D/21}\mathrm{\Gamma }(D/2)}.$$
(15)   
Therefore, if the probability density $`f(\alpha )`$ do not allow $`\alpha `$ to be smaller than $`D`$, the only imaginary term in Eq.(14) comes from $`a(E+iϵ)`$, and it is of order $`p^2E^{3/2}`$, yielding a value of order $`p^5`$ at the peak. On the other hand, if $`\alpha `$ can be smaller than $`D`$, an imaginary part of order $`p^2f(D)`$ arises from the pole.

Led by the functional form of the self-energy in Eq.(14), one can consider self-energies of the form $`f(E)E_0(p)`$, also when $`\rho `$ is not small, and the single defect approximation no longer holds. This is the idea lying under the well-known CPA approximation , where one sets $`\mathrm{\Sigma }(z,p)=(\mathrm{\Gamma }(z)1)E_0(p)`$. A self consistency equation can be readily written:

$$\overline{\frac{1+\alpha \mathrm{\Gamma }(z)}{\mathrm{\Gamma }(z)a\left(z/\mathrm{\Gamma }(z)\right)\left(1+\alpha \mathrm{\Gamma }(z)\right)}}=0.$$
(16)   
It is clear that the width of the $`S(q,E)`$ critically depends on the value of $`\mathrm{\Gamma }(0+iϵ)`$. For the flat distribution of $`\alpha `$ introduced in Eq.(5), one can solve Eq.(16) in the limit of small $`\rho `$, as $`\mathrm{\Gamma }(0)=1+b\rho +𝒪(\rho ^2)`$ obtaining

$$b=D\left(1D\mathrm{log}\frac{D}{|\lambda +D|}+i\pi \frac{D}{\lambda }\theta (D+\lambda )\right),$$
(17)   
where $`\theta (x)=1`$ for $`x<0`$ and zero otherwise. It is easy to check that the single-defect result is exactly the same as in Eq.(17).

Fixing from now on $`D=3`$, Eq.(16) in the general case, can be numerically solved. We choose to write it as a fixed-point equation, and solve it recursively. The only tricky part, is the evaluation of $`a(z)`$ defined in Eq.(12), that we make by calculating the unperturbed density of states by a Monte Carlo simulation. In this way, we are able to obtain estimates of $`\mathrm{\Gamma }(z)`$ with an accuracy of $`10^4`$. We find that in the CPA approximation, there are the same two regimes as in the single-defect computation, separated by a critical line, whose coordinates are shown in table I. The appearance of the $`p^5`$ regime coincides, and it is due to, the vanishing of the negative-energy spectrum. However, as explained in the introduction, one expects rather that what vanish is the order $`\rho `$ contribution to the density of states in this region. In fact, we expect a non-vanishing density of states all the way down to $`\lambda =1`$, where no negative spring exists. It is remarkable that when $`\rho =1`$, the CPA approximation keeps a non-vanishing fraction of negative eigenvalues up to $`\lambda =1.2`$, quite close to the correct value.

The agreement between the CPA approximation and the results from numerical simulations for the $`S(E,p)`$ turns out to be better than $`5\%`$ in the two extreme cases (see figs. 2,3). That is the one where there is no negative springs ($`\lambda =1`$) and the one where the spring constant can be very negative ($`\lambda =10`$). Close to the CPA critical line (in the unstable crystal side) the agreement is still quite good. On the (CPA) stable border side, a scaling between $`p^6`$ and $`p^5`$ is found, depending on the density of defects, $`\rho `$. For example, with the value of the density of defects $`\rho =0.1`$, in fig. 1) the transition between the two regimes found by the CPA approximation at $`\lambda =2.15`$ is very evident.

We numerically computed the structure factor $`S(E,q)`$ for our model utilizing the method of moments , which allows to study the statistical properties of the eigenvalues and the eigenvectors of large dynamical matrices, avoiding their diagonalization. This method is a clever modification of the Lanczos method and it shares the same weakness, namely the lack of orthogonality when a too large number of moments is computed. Another limitation is the necessity of setting a finite value of $`ϵ`$ in Eq.(4). A reasonable value for $`ϵ`$ is the mean distance, between eigenstates, which is roughly given by $`(2D\lambda )/V`$, where $`V`$ is the lattice volume (that is $`10^5`$ in our $`96^3`$ lattice). If the width of the peak is comparable with $`ϵ`$, the results will be definitely affected by finite-size effects. The limitation related with the number of moments, is not serious for the central part of the $`S(q,E)`$ curve, but can be rather strong if one wants to calculate the tails of the distribution. In practice, we have used $`30`$ moments for $`10`$ different disorder realizations, finding very satisfying results, unless the peak height approaches values of order $`10^5`$, when finite size effects turns out to be important. We fit our results to the Breit-Wigner form:

$$S(E,q)=\frac{𝒩}{(EE_0)^2+\mathrm{\Sigma }^2},$$
(18)   
which satisfactorily describes the peak in all cases, although usually overestimates the tails of the distribution. The position of the peak $`E_0`$ is linear in $`p^2`$ as expected. More importantly, the results from numerical simulations (fig. 4) show very clearly that the critical line previously discussed is an artifact of the CPA approximation and a real system actually becomes unstable when there is an extensive number of negative springs, no matter how small is the fraction of them, as expected from the simple analytic considerations sketched previously. For larger values of $`p`$ the $`p^5`$ contribution is no longer negligible and a complicated intermediate scaling appears.

In this work, we used the fact that both the single-defect approximation and the CPA approximation  allow to write the self-energy of a disordered harmonic solid, for small $`E`$ and $`p`$, as $`f(E)E_0(p)`$. An unavoidable consequence of this functional form is that violations of the Rayleigh $`p^5`$ scaling of the width of the peaks of the $`S(p,E)`$ are intrinsically tied with the appearance of a negative energy spectrum (i. e. solid instability). We have argued that an extensive number of negative eigenvalues should appear, as soon as the spring-constants are allowed to be negative. Our numerical calculations on $`96^3`$ disordered lattices confirm the above picture. For wave-lengths in the range $`[a,10a]`$ ($`a`$ being the lattice constant), a crossover regime may come out, depending on the strength of the disorder, from the competition between the $`p^2`$ term and $`p^5`$ term in the self-energy. Therefore, the vibrational excitations of the disordered solid can look qualitatively similar to the finite temperature Instantaneous Normal Modes of supercooled liquids and glasses, where negative eigenvalues are always present.

We gratefully acknowledge interesting discussions with A. Gonzalez, M. Mézard, G. Ruocco and G. Viliani. V.M.M. is a M.E.C. fellow and has been partially suported by CYCyT(AEN97-1708 and AEN99-1693). Our numerical computations have been carried out on the Kalix2 pentium cluster of the University of Cagliari.