format: text

text:
# The dynamical structure factor in topologically disordered systems

## I Introduction

Inelastic X-ray scattering (IXS) experiments  and inelastic neutron scattering  on structural glasses and supercooled liquids provided useful information on the dynamics of their amorphous structure, at frequencies larger than $`0.1`$ THz. Those experiments show a regime, when the wavelength of the plane wave is comparable with the inter-particle distance, where the vibrational spectrum can be understood in terms of propagation of quasi-elastic sound waves, the so-called high frequency sound. This high-frequency sound has also been observed in molecular dynamical simulations of strong  and fragile  liquids, and it displays several rather universal features  (see however Ref.  for a dissenting view regarding the silica properties). In particular, a peak is observed in the dynamical structure factor at a frequency which depends linerly on the exchanged momentum $`p`$, in the region $`0.1p_01.0p_0`$, $`p_0`$ being the position of the first maximum in the static structure factor. When extrapolated to zero momentum, this linear dispersion relation yields the macroscopic speed of sound. The width of the spectral line, $`\mathrm{\Gamma }`$ is well fitted by

$$\mathrm{\Gamma }(p)=Ap^x,x2,$$
(1)
with $`A`$ displaying a very mild (if any) temperature dependence. Although the same scaling of the spectral line is found in hydrodynamics , in this computation the proportionality constant is basically the viscosity, which shows a very strong temperature dependence which is not observed. Moreover, the same scaling of $`\mathrm{\Gamma }`$ has been found in harmonic Lennard-Jones glasses , and one can safely conclude that the $`p^2`$ broadening of the high-frequency sound has a physical origin different from hydrodynamics. Other interesting problems related with the high-frequency vibrational excitations of these topologically disordered systems  regard the origin of the Boson peak or the importance of localization properties to understand the dynamics of supercooled liquids .

The variety of materials in which the $`p^2`$ broadening appears suggests a straightforward physical motivation. However, the simplest conceivable approximation, a wave propagating on an elastic medium in the presence of random scatterers, yields Rayleigh dispersion: $`\mathrm{\Gamma }p^4`$. This result is very robust: as soon as one assumes the presence of an underlying medium on which the sound waves would propagate undisturbed, as in the disordered-solid model , the $`p^4`$ scaling appears even if one studies the interaction with the scatterers non-perturbatively . In this paper we want to show that when the distinction between the propagating medium and the scatterers is meaningless (as it happens for topologically disordered systems), the $`p^2`$ scaling is recovered. Recently, it has been shown  that the numerical solution of the mode coupling equations, modified for the study of the glassy phase, describes in a qualitatively correct way the range of frequencies explored by IXS scattering, including the high frequency sound and the boson peak. In that theoretical approach though, one cannot find a clear scaling law, $`\mathrm{\Gamma }`$ being proportional to $`p^2`$ only at very low momentum and following a more complicated law for higher values.

We want to investigate the problem from the point of view of statistical mechanics of random matrices , by assuming that vibrations are the only motions allowed in the system. The formalism we shall introduce, however, is not limited to the investigation of the high frequency sound and it could be straightforwardly applied in different physical contexts.

Let us look more carefully at the relation between the high frequency behaviour and vibrational dynamics in glasses and at the relation between vibrations and random matrices. The dynamical structure factor for a system of $`N`$ identical particles is defined as :

$$S(p,\omega )=\frac{1}{N}\underset{i,j}{}𝑑t\mathrm{e}^{\mathrm{i}\omega t}\mathrm{e}^{\mathrm{i}p\left(r_j(t)r_i(0)\right)},$$
(2)
where $`\mathrm{}`$ denotes the average over the particles positions $`r_j`$ with the canonical ensemble.

A well known approach to the calculation of the $`S(p,\omega )`$ in the high-frequency (THz) region, consists in taking into account only the vibrational modes of the system. The basic assumption of this calculational strategy is that there is a clear separation among the time scales of the ’fast’ degrees of freedom, which are supposedly modeled as vibrations, and the ’slow’ ones, which give rise to the diffusion of the particles. One expects this assumption to be reasonably accurate in the glass and in the super-cooled liquid phase.

When focusing only on the short time region, the slow motions can be thought as ’frozen’, hence one studies the displacements $`u`$ around ’quenched’ positions $`x`$, by writing the position of the $`i`$-th particle as $`r_i(t)=x_i+u_i(t)`$, and linearizing the equations of motion. Then one is naturally lead to consider the spectrum of the Hessian matrix of the potential, evaluated on the positions $`x`$ of the particles. Calling $`\omega _n^2`$ the eigenvalues of the Hessian matrix and $`e_n(i)`$ the corresponding eigenvectors, the one excitation approximation to the $`S(p,\omega )`$ at non zero frequency is given in the classical limit by:

$`S^{(1)}(p,\omega )`$ $`=`$ $`{\displaystyle \frac{k_\mathrm{B}T}{m\omega ^2}}{\displaystyle \underset{n=1}{\overset{N}{}}}Q_n(p)\delta (\omega \omega _n),`$ (3)
$`Q_n(p)`$ $`=`$ $`|{\displaystyle \underset{i}{}}pe_n(i)\mathrm{e}^{\mathrm{i}px_i}|^2.`$ (4)
However, in the supercooled liquid phase, one cannot always assume that in a typical configuration at equilibrium all the normal modes have positive eigenvalues, negative eigenvalues representing the situation where some particles are moving away from the position $`x`$. This observation has motivated the so-called instantaneous normal modes (INM) approach  in the study of supercooled liquids. In that approach the normal modes in a typical configuration are supposed to describe the short time dynamics of the particles in the liquid phase, while it has been suggested  that diffusion properties could be studied by considering the localization properties of the normal modes of negative eigenvalues.

On the other hand, it has been argued that in a very broad class of glassy systems below a critical temperature , called dynamical temperarture $`T_D`$ (for spin glasses) or mode coupling temperature $`T_{MCT}`$ (for structural glasses), the particles spend the most of their time around the minima of the potential energy surface. Recent numerical simulations seem to confirm that picture .

As a consequence, in the glassy phase the high frequency window ($`.1`$ Thz) of the dynamical structure factor is supposed to be well described by the propagation of quasi-harmonic excitations around a quenched structure . By assuming that this structure corresponds to one of the many minima of the potential energy, one can introduce a harmonic approximation where only the vibrations around these minima are considered, and all the dynamical information is encoded in the spectral properties of the Hessian matrix evalutated on the rest positions, where it has not negative eigenvalues. It has been shown using molecular dynamics simulation that below the experimental glass transition temperature the thermodynamical properties of typical strong glasses  are in a good agreement with such an assumption. Let us note that, within this assumption, the dynamical structure factor (4) depend on the minimum around which the particles oscillate, hence the knowledge of an infinite number of disordered positions $`x`$ seems to be needed. Here we shall make the reasonable assumption that the dynamical structure factor is self-averaging, a type of assumption which generally turns out to be correct for macroscopic observables in disordered systems.

Hence we have turned our problem into the study of spectral properties of a kind of random matrix, once the probability distribution of the elements of the matrix is given. More precisely, the problem of the high-frequency dynamics of the system can be reduced, in its simplest version, to the consideration of a very particular type of random matrices, the euclidean matrices , in which the entries are deterministic functions (the derivatives of the potential) of the random positions of the particles. If the system has a conservation law such as momentum conservation in our case, due to translational invariance, all the realizations of the random matrix have a common normal mode with zero eigenvalue: the uniform translation of the system. An euclidean matrix is determined by the particular deterministic function that we are considering, and by the probabilistic distribution function of the particles. However, for the sake of simplicity we shall concentrate here on the simplest kind of euclidean matrices . We consider $`N`$ particles placed at positions $`x_i,i=1,2,\mathrm{},N`$ inside a box of volume $`V`$, where periodic boundary conditions are applied. Then, the scalar euclidean matrices are

$$M_{ij}\delta _{ij}\underset{k=1}{\overset{N}{}}f(x_ix_k)f(x_ix_j),i,j=1,2\mathrm{},N.$$
(5)
In the above expression, $`f(x)`$ is an arbitrary but deterministic scalar function depending on the distance between pairs of particles, and the positions $`\{x\}`$ of the particles are drawn by a probabilistic law $`P[x]`$. Notice that the matrix (5) preserves translation invariance, since the uniform vector $`e_0(i)=const`$ is an eigenvector of $`M`$ with zero eigenvalue. Let us note furthermore that there are not internal indices, since for the sake of simplicity the particle displacements are restricted to be all collinear. Of course, in the study of the vibrations of glasses, one should take into account the vectorial nature of those vibrations and split the Hessian matrix in its transversal and longitudinal (with respect to the directions of the external momentum $`\stackrel{}{p}`$) parts. A careful analysis of that feature will be presented elsewhere .

Simple as these matrices are, many physical problems ranging from vibrational spectra in amorphous systems and electron hopping in amorphous semiconductors , to instantaneous normal modes (INM) in liquids  and combinatorial optimization , can be described with them using an appropriate choice of $`P`$. In this paper we will study the simplest case, in which $`P`$ is a uniform distribution function.

The dynamical structure factor for a scalar euclidean-matrix is given by

$`S_E(p,E)`$ $`=`$ $`\overline{{\displaystyle \underset{n}{}}Q_n(p)\delta (EE_n)},`$ (6)
$`Q_n(p)`$ $`=`$ $`{\displaystyle \frac{1}{N}}\left|{\displaystyle \underset{i=1}{\overset{N}{}}}e_n(i)\mathrm{e}^{\mathrm{i}px_i}\right|^2,`$ (7)
$`S(p,\omega )`$ $`=`$ $`2\omega S_E(p,\omega ^2),`$ (8)
where the overline stands for the average over the particles position and we have given the definition either in the eigenvalue space ($`S_E(p,E)`$) and in the frequency space ($`S(p,\omega )`$).

Finally, let us mention that not all the topologically disordered systems are represented by an euclidean random matrix. For instance, diffusion on random graphs  and on small-world network  where the connecting sites are chosen with a probabilistic law, are described by the Laplacian operator $`𝒲_{ij}`$. The matrix element $`𝒲_{ij}`$, $`ij`$ is one if $`i`$ and $`j`$ are connected and $`0`$ otherwise, while the diagonal part is tailored to have $`_i𝒲_{ij}=0`$, as required by the conservation of the number of diffused particles.

The layout of the rest of this paper is organized as follows. In section II we shall apply an expansion in the inverse of the particle density formulation to the problem of calculating the $`S_E(p,E)`$. The same calculation can be performed using the field theory of Ref. , adapted to this problem , with identical results . The analytical computations will be confronted with a numerical simulation in section III. In section IV we shall briefly consider the problem of including particle correlations on our calculation, which shall produce some qualitatively new results, regarding the exponential tail of the density of states. In section V we shall present our conclusions.

## II The 1/$`\rho `$ expansion

The basic object that we will calculate is the resolvent

$`G(p,z)`$ $`=`$ $`{\displaystyle \frac{1}{N}}\overline{{\displaystyle \underset{ij}{}}e^{ip(x_ix_j)}\left[{\displaystyle \frac{1}{zM}}\right]_{ij}},`$ (9)
from which the dynamical structure factor is recovered using the distribution identity $`1/(x+i0^+)=P(1/x)\mathrm{i}\pi \delta (x)`$:

$$S_E(p,E)=\frac{1}{\pi }\underset{\eta 0^+}{lim}\mathrm{I}mG(p,E+\mathrm{i}\eta ).$$
(10)
The resolvent will also give us access to the density of states of the system:

$`g_E(E)`$ $`=`$ $`{\displaystyle \frac{1}{N\pi }}\underset{\eta 0^+}{lim}\mathrm{I}m\overline{{\displaystyle \underset{i=1}{\overset{N}{}}}\left[{\displaystyle \frac{1}{E+\mathrm{i}\eta M}}\right]_{ii}},`$ (11)
$`=`$ $`{\displaystyle \frac{1}{\pi }}\underset{\eta 0^+}{lim}\mathrm{I}m\underset{p\mathrm{}}{lim}G(p,E+\mathrm{i}\eta ).`$ (12)
We shall set up a perturbative expansion in the inverse density of particles, $`1/\rho `$. The zeroth order of this expansion (to be calculated in subsection II A) corresponds to the limit of infinite density, in which the system is equivalent to an elastic medium. In this limit the resolvent turns out to be extremely simple:

$`G(p,z)`$ $`=`$ $`{\displaystyle \frac{1}{zϵ(p)}},`$ (13)
$`ϵ(p)`$ $`=`$ $`\rho [\stackrel{~}{f}(0)\stackrel{~}{f}(p)].`$ (14)
In the above expression $`\stackrel{~}{f}(p)`$ is the Fourier transform of the function $`f`$ defined in Eq.(5), that due to its spherical symmetry, behaves like $`ϵ(p)=𝒪(p^2)`$ for small $`p`$. We see that the dynamical structure function has two delta functions at frequencies linear in $`p`$: $`\omega =\pm \sqrt{ϵ(p)}`$. It is then clear that Eq.(13) represents the undamped propagation of sound in our elastic medium, with a dispersion relation controlled by the function $`ϵ(p)`$ . The order $`1/\rho `$ corrections to Eq.(13) will be calculated in subsection II B. They take the form of a complex self-energy, $`\mathrm{\Sigma }(p,z)`$, yielding

$`G(p,z)`$ $`=`$ $`{\displaystyle \frac{1}{zϵ(p)\mathrm{\Sigma }(p,z)}},`$ (15)
$`S_E(p,E)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle \frac{Im\mathrm{\Sigma }(p,E)}{\left(Eϵ(p)Re\mathrm{\Sigma }(p,E)\right)^2+\left(Im\mathrm{\Sigma }(p,E)\right)^2}}.`$ (16)
The dynamical structure factor, is no longer a delta function, but close to its maxima it has a Lorentzian shape. From Eq.(16) we see that the real part of the self-energy renormalizes the dispersion relation, and thus it should behaves proportionallyto $`p^2`$ at small momentum. The width of the spectral line is instead controlled by the imaginary part of the self-energy, and at this order in $`1/\rho `$, it scales like $`p^4`$ for small momentum: at this level of accuracy we recover Rayleigh scattering. The analytic calculation to this order can be confronted with numerical simulations (see section III), and it turns-out to be rather accurate for momenta of the order of the inverse interparticle distance. The density of states for small values of $`\omega `$ is also well reproduced. The maybe most interesting result of this paper is obtained in section II C, where the calculation of the self-energy is pushed up to order $`1/\rho ^2`$. The order $`1/\rho ^2`$ contribution to the spectral line-width is proportional to $`p^2`$ for small momentum, which is a qualitatively new feature. Even if the density of the system is high, for small-enough momenta the two loops contribution to the imaginary part of the self-energy is dominating. We believe this to be the underlying physical reason for the widespread experimental finding of $`p^2`$ scaling of the width of the dynamical structure factor. On the other hand, for larger momenta the order $`1/\rho `$ contribution is dominating, but it grows with $`p`$ significantly more slowly than $`p^4`$. This maybe of some importance in the experimentally relevant momentum range.

### A The leading order

Our strategy for computing the resolvent (9) is straightforward: we consider it as the sum of a geometric series. One expands in series of $`1/z`$ the Fourier transform of the resolvent, see Eq.(9), in the following way:

$$G(p,z)=\frac{1}{z}\underset{R}{}\left(\frac{1}{z}\right)^R\overline{M^R(p)}$$
(17)
with

$`M^R(p)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{k_0,k_1\mathrm{}k_R}{}}e^{ipx_{k_0}}\left(\delta _{k_0,k_1}{\displaystyle \underset{z_1}{}}f(x_{k_0}x_{z_1})f(x_{k_0}x_{k_1})\right)\mathrm{}`$ (18)
$`\mathrm{}`$ $`\left(\delta _{k_{R1},k_R}{\displaystyle \underset{z_R}{}}f(x_{k_{R1}}x_{z_N})f(x_{k_{R1}}x_{k_R})\right)e^{ipx_{k_R}}`$ (19)
Since we are interested in the singularities (branch-cuts or poles) on the complex $`z`$ plane of the resolvent $`G(p,z)`$, see Eqs.(13,15), it will be crucial to sum the series in Eq.(19) at all orders in $`1/z`$. From now on we shall assume that the function $`f`$ has a Fourier transform (denoted by $`\stackrel{~}{f}`$ in the following) which is a rather severe limitation. In section IV we will see that, provided that particle correlations are taken into account on the calculation, this limitation should not bother us. One can readily check from Eq.(19) that the contribution of the terms where the two particle labels in the argument of any $`f`$ function are the same, vanish. This can be taken care of automatically by setting $`f(0)=0`$. Since we have changed the value of the function at a single point, this will have no consequences while averaging over the positions of different particles.

The very first step is to split the parenthesis in Eq.(19), into $`2^R`$ sums (we keep the Kronecker deltas as a reminder of which $`k_i`$ index should be equal to which):

$`M^R(p)`$ $`=`$ $`{\displaystyle \underset{k_0,z_0,\mathrm{},z_{R1}}{}}{\displaystyle \frac{e^{ip(x_{k_0}x_{k_R})}}{N}}[\delta _{k_0,k_1}f(x_{k_0}x_{z_0})][\delta _{k_1,k_2}f(x_{k_1}x_{z_1})]\mathrm{}`$ (21)
$`\mathrm{}[\delta _{k_{R1},k_R}f(x_{k_{R1}}x_{z_{R1}})]+`$
$`+`$ $`{\displaystyle \underset{k_0,k_1,z_1,\mathrm{},z_{R1}}{}}{\displaystyle \frac{e^{ip(x_{k_0}x_{k_R})}}{N}}[f(x_{k_0}x_{k_1})][\delta _{k_1,k_2}f(x_{k_1}x_{z_1})]\mathrm{}`$ (24)
$`\mathrm{}[\delta _{k_{R1},k_R}f(x_{k_{R1}}x_{z_{R1}})]+`$
$`\mathrm{}`$
$`+`$ $`{\displaystyle \underset{k_0,k_1,\mathrm{},k_R}{}}{\displaystyle \frac{e^{ip(x_{k_0}x_{k_R})}}{N}}[f(x_{k_0}x_{k_1})][f(x_{k_1}x_{k_2})]\mathrm{}[f(x_{k_{R1}}x_{z_{R1}})]`$ (25)
The disorder average in Eq.(25) is simply obtained by performing the integrals:

$$\overline{M^R(p)}=\underset{j}{}\frac{d^Dx_j}{V}M^R(p)$$
(26)
Let us estimate the order of magnitude of the generic term in Eq.(25), after the disorder average is performed. If there are no index repetitions, we will have an overall $`V^R`$ factor from the disorder average, because although there are $`R+1`$ independent particles the translational invariance allows one to eliminate one integral. We have $`N!/(NR1)!N^{R+1}`$ ways of choosing the particles labelling, and an extra $`N^1`$ factor from the definition of the resolvent (9). Therefore, each of the above $`2^R`$ sums will be proportional to $`\rho ^R`$ ($`\rho =N/V`$ is the particle density). When two particle label coincide, we will have a missing $`N`$ factor from the sums, and another missing $`V^1`$ factor from the position average, and the sum will be of order $`\rho ^{R1}`$. It is thus clear that the $`\overline{M^R(p)}`$ will be a polynomial of order $`R`$ in $`\rho `$, with no term of order $`\rho ^0`$. In fact, the lowest order in $`\rho `$ arise when there are only two different particle labels on the sums of Eq.(25), which is easily seen to be linear in $`\rho `$:

$$\overline{M^R(p)}=\rho ^R_R^{(R)}(p)+\rho ^{R1}_{R1}^{(R)}(p)+\mathrm{}+\rho _1^{(R)}(p).$$
(27)
Since $`M^R(p=0)=0`$, for every particle configuration, all the coefficient $`_k^{(R)}(p)`$ in Eq.(27) vanish at zero external momentum, $`p=0`$.

Our zeroth-order approximation will consist in keeping only the leading term

$$\overline{M^R(p)}\rho ^R_R^{(R)}(p),$$
(28)
the first $`1/\rho `$ correction will consist in setting

$$\overline{M^R(p)}\rho ^R_R^{(R)}(p)+\rho ^{R1}_{R1}^{(R)}(p),$$
(29)
and so on.

Let us now calculate $`_R^{(R)}(p)`$. As for any multidimensional integral, choosing with care the order in which the integrations are performed can tremendously simplify the task. In the case of Eq.(25) the matrix product yields a fairly natural ordering: from left to right, and when we encounter a diagonal term ($`\delta _{k_i,k_{i+1}}`$) we shall first average the position of the $`z_i`$ particle. One easily shows that (recall that the Kronecker delta is simply a reminder, since $`k_i`$ is actually equal to $`k_{i+1}`$):

$`{\displaystyle 𝑑x_{z_i}e^{ipx_{k_i}}\delta _{k_i,k_{i+1}}f(x_{k_i}x_{z_i})}`$ $`=`$ $`e^{ipx_{k_{i+1}}}\stackrel{~}{f}(0),`$ (30)
$`{\displaystyle 𝑑x_{k_i}e^{ipx_{k_i}}[f(x_{k_i}x_{k_{i+1}})]}`$ $`=`$ $`e^{ipx_{k_{i+1}}}[\stackrel{~}{f}(p)].`$ (31)
Therefore every integral simply shifts the external moment, $`p`$, one step forward, and leaves behind a $`\stackrel{~}{f}(0)`$ or a $`\stackrel{~}{f}(p)`$. When one arrives to the $`R+1`$ particle, there is no integral left to be done, and both exponentials cancel. Then we find that

$$\rho ^R_R^{(R)}(p)=\left[\rho \stackrel{~}{f}(0)\rho \stackrel{~}{f}(p)\right]^R,$$
(32)
and therefore we obtain the result already anticipated in Eq.(13)

$$G_0(p,z)=\frac{1}{zϵ(p)}.$$
(33)
We see that at leading order, a plane wave with momentum $`p`$ is actually an eigenstate of the matrix $`M`$ with eigenvalue $`E`$, and the disorder does not play any relevant role. In other words, inside a wavelength $`2\pi /p`$ there is always an infinite number of particles, smoothing out the density fluctuations of the particles: the system reacts as an elastic medium.

Let us finally obtain the density of states at this level of accuracy, using Eq(12) and $`G_0`$:

$$g_E(E)=\delta \left(E\rho \stackrel{~}{f}(0)\right).$$
(34)
We obtain a single delta function at $`\rho \stackrel{~}{f}(0)`$, which is somehow contradictory with our result for the dynamical structure factor: from the density of states one would say that the dispersion relation is Einstein’s like, without any momentum dependence! The way out of this contradiction is of course that in the limit of infinite $`\rho `$ both $`ϵ(p)`$ and $`\rho \stackrel{~}{f}(0)`$ diverge. The delta function in eq. (34) is the leading term in $`\rho `$, while the states which contribute to the dynamical structure factor appear only in the subleading terms in the density of states.

### B One loop

For the calculation of the $`_{R1}^{(R)}(p)`$ term, we need to consider only one particle label repetition. In other words, if we call $`N_R`$ the average number of particles within the interaction range of the potential, we are calculating the 1/$`N_R`$ correction of the formula $`\overline{M^R}=\overline{M}^R`$, which neglects the statistical fluctuation of the matrix $`M`$ due to disorder.

The easiest way to proceed is as follows: one first identifies the places where the particle label repetition arises in the sequence of particle labels

$$(\mathrm{}1\mathrm{}1\mathrm{}),$$

and then one shifts the moment using Eqs.(30) until the place of the first repeated particle label. Then, one applies the backward version of the moment shift idea, from $`k_R`$ until the second particle label repetition, and one finally integrates over the position of the repeated particle. In this way one gets:

$`{\displaystyle \frac{\rho ^{R1}_{R1}^{(R)}(p)}{z^{R+1}}}`$ $`=`$ $`{\displaystyle \underset{a+b+c=R2}{}}{\displaystyle \frac{[\rho (\stackrel{~}{f}(0)\stackrel{~}{f}(p))]^a}{z^{a+1}}}{\displaystyle \frac{[\rho (\stackrel{~}{f}(0)\stackrel{~}{f}(p))]^c}{z^{c+1}}}\times `$ (35)
$`\times `$ $`{\displaystyle \frac{dq}{\rho (2\pi )^D}[\rho (\stackrel{~}{f}(pq)\stackrel{~}{f}(q))]^2\frac{[\rho (\stackrel{~}{f}(0)\stackrel{~}{f}(q))]^b}{z^{b+1}}},`$ (36)
which is easily seen to yield

$$\underset{R=0}{\overset{\mathrm{}}{}}\frac{\rho ^{R1}_{R1}^{(R)}(p)}{z^{R+1}}=[G_0(z,p)]^2\frac{1}{\rho }\frac{dq}{(2\pi )^D}G_0(q)[\rho (\stackrel{~}{f}(pq)\stackrel{~}{f}(q))]^2.$$
(37)
We want to obtain the self-energy term in (15), which will receive contributions of all powers in $`1/\rho `$:

$$\mathrm{\Sigma }(p,z)=\mathrm{\Sigma }_1(p,z)+\mathrm{\Sigma }_2(p,z)+\mathrm{}$$
(38)
The relation between the above expresion and the resolvent is given by the Dyson resummation:

$$G(p,z)G_0(p,z)=G_0(p,z)\left[\mathrm{\Sigma }_1(p,z)+G_0(p,z)\mathrm{\Sigma }_1^2(p,z)+\mathrm{\Sigma }_2(p,z)+\mathrm{}\right]G_0(p,z).$$
(39)
Comparing Eqs.(39) and (37), we conclude that

$$\mathrm{\Sigma }_1(p,z)\frac{1}{\rho }\frac{d^dq}{(2\pi )^d}G_0(q,z)\left[\rho \stackrel{~}{f}(p\text{q})\rho \stackrel{~}{f}(\text{q})\right]^2.$$
(40)
The $`1/\rho ^2`$ contributions to the self-energy will provide the $`G_0^3(p,z)\mathrm{\Sigma }_1^2(p,z)`$ term of Eq.(39), and also the more physically interesting $`\mathrm{\Sigma }_2(p,z)`$ term.

Let us study in details the low exchanged momentum limit of Eq.(40). It is clear that at $`p=0`$ the self-energy vanishes, as required by the translational invariance. We need to expand $`\stackrel{~}{f}(pq)`$ for small $`p`$, which due to the spherical symmetry of $`\stackrel{~}{f}`$ yields

$`\stackrel{~}{f}(pq)`$ $`=`$ $`\stackrel{~}{f}(q)(pq){\displaystyle \frac{\stackrel{~}{f}^{}(q)}{q}}+𝒪(p^2),`$ (41)
$`=`$ $`\stackrel{~}{f}(q)+(pq){\displaystyle \frac{ϵ^{}(q)}{q\rho }}+𝒪(p^2),`$ (42)
where we have used Eq.(14) for the derivative of $`\stackrel{~}{f}`$. Substituting (42) in (40), and performing explicitly the trivial angular integrations in dimensions $`d`$ we obtain

$`\mathrm{\Sigma }_1(p,z)`$ $``$ $`p^2{\displaystyle \frac{2^{1d}}{\rho d\pi ^{d/2}\mathrm{\Gamma }(d/2)}}{\displaystyle _0^{\mathrm{}}}𝑑qq^{d1}{\displaystyle \frac{\left[ϵ^{}(q)\right]^2}{zϵ(q)}},`$ (43)
$`=`$ $`p^2{\displaystyle \frac{2^{1d}}{\rho d\pi ^{d/2}\mathrm{\Gamma }(d/2)}}{\displaystyle _0^{ϵ(q=\mathrm{})}}𝑑ϵ{\displaystyle \frac{\left[q(ϵ)\right]^{d1}}{q^{}(ϵ)(zϵ)}}.`$ (44)
In the last equation, we have denoted with $`q(ϵ)`$ the inverse of the function $`ϵ(q)`$. Setting now $`z=E+\mathrm{i0}^+`$, and observing that $`ϵ(p)Ap^2`$ for small $`p`$, we readily obtain

$`Re\mathrm{\Sigma }_1(p,E+\mathrm{i0}^+)`$ $``$ $`p^2{\displaystyle \frac{2^{1d}}{\rho d\pi ^{d/2}\mathrm{\Gamma }(d/2)}}P{\displaystyle _0^{ϵ(q=\mathrm{})}}𝑑ϵ{\displaystyle \frac{\left[q(ϵ)\right]^{d1}}{q^{}(ϵ)(Eϵ)}},`$ (46)
$`Im\mathrm{\Sigma }_1(p,E+\mathrm{i0}^+)`$ $``$ $`{\displaystyle \frac{\pi 2^{2d}A}{\rho d\pi ^{d/2}\mathrm{\Gamma }(d/2)}}p^2[q(E)]^d.`$ (47)
Since the principal part is a number of order one, the real part of the self-energy scales like $`p^2`$ (possibly with logarithmic corrections), and thus the speed of sound of the system renormalizes due to the $`1/\rho `$ corrections. As a consequence, the function $`q(E)`$ is proportional to $`E^{1/2}p`$ at the maximum of the function of $`p`$ $`S_E(p,E)`$, and the width of the peak of the $`S_E(p,E)`$ will scale like $`p^{d+2}`$. It is then easy to check (see (8)) that in frequency space the width of the spectral line will scale like

$$\mathrm{\Gamma }p^{d+1},$$
(48)
as one would expect from Rayleigh scattering considerations. The result (48) for the asymptotic regime $`p<<1`$ has been found at the one loop level. In order to predict correctly the spectral properties at very low external momentum $`p`$, it turns out that one must study the behaviour of the two loop contribution, as we shall see in the next section. Nevertheless, the one loop result is already a good starting point to perform detailed comparisons with the numerical simulations. To this end it will be useful in the following to introduce a concrete example with a specific choice of the function $`f(x)`$, allowing a complete analytical computation.

#### 1 The Gaussian case

Our choice for the function $`f`$ shall be the following

$`f(x)`$ $`=`$ $`\mathrm{e}^{\frac{x^2}{2\sigma ^2}},`$ (49)
$`\stackrel{~}{f}(p)`$ $`=`$ $`(2\pi \sigma ^2)^{d/2}\mathrm{e}^{\sigma ^2\frac{p^2}{2}}.`$ (50)
The parameter $`\sigma `$ sets our length-scale, and for phenomenological purposes might be identified with the inverse of the first maximum of the static structure factor. With this choice, we have

$`ϵ_{max}`$ $``$ $`\underset{q\mathrm{}}{lim}ϵ(q)=\rho (2\pi \sigma ^2)^{d/2},`$ (51)
$`ϵ(q)`$ $`=`$ $`ϵ_{max}\left(1\mathrm{e}^{\sigma ^2\frac{q^2}{2}}\right),`$ (52)
$`q(ϵ)`$ $`=`$ $`{\displaystyle \frac{1}{\sigma }}\sqrt{2\mathrm{ln}{\displaystyle \frac{ϵ_{max}}{ϵ_{max}ϵ}}}.`$ (53)
It is then straightforward to obtain in three dimensions

$`\mathrm{\Sigma }_1(p,z)`$ $`=`$ $`\sqrt{{\displaystyle \frac{2\sigma ^2}{\pi }}}{\displaystyle _0^{ϵ_{max}}}𝑑ϵ{\displaystyle \frac{𝒢(q(ϵ),p)}{zϵ}},`$ (54)
$`𝒢(q,p)`$ $`=`$ $`q\mathrm{e}^{\sigma ^2\frac{q^2}{2}}\left[1+\mathrm{e}^{\sigma ^2p^2}{\displaystyle \frac{\mathrm{sh}(2\sigma ^2pq)}{2\sigma ^2pq}}2\mathrm{e}^{\sigma ^2\frac{p^2}{2}}{\displaystyle \frac{\mathrm{sh}(\sigma ^2pq)}{\sigma ^2pq}}\right]`$ (55)
Therefore, setting $`z=E+\mathrm{i0}^+`$, we find

$`Re\mathrm{\Sigma }_1(p,E+\mathrm{i0}^+)`$ $`=`$ $`\sqrt{{\displaystyle \frac{2\sigma ^2}{\pi }}}P{\displaystyle _0^{ϵ_{max}}}𝑑ϵ{\displaystyle \frac{𝒢(q(ϵ),p)}{Eϵ}},`$ (56)
$`Im\mathrm{\Sigma }_1(p,E+\mathrm{i0}^+)`$ $`=`$ $`\sqrt{2\pi }\sigma ^2𝒢(q(E),p).`$ (57)
One can easily check that the imaginary part of the self-energy on the peak is of order $`p^{d+2}`$, as previously announced.

Turning now to the density of states, using Eq.(12), Eqs.(56) and (57), simplify now to

$`Re\mathrm{\Sigma }_1(\mathrm{},E+\mathrm{i0}^+)`$ $`=`$ $`\sqrt{{\displaystyle \frac{2\sigma ^2}{\pi }}}P{\displaystyle _0^{ϵ_{max}}}𝑑ϵ{\displaystyle \frac{q(ϵ)\mathrm{e}^{\sigma ^2q(ϵ)^2/2}}{Eϵ}},`$ (58)
$`Im\mathrm{\Sigma }_1(\mathrm{},E+\mathrm{i0}^+)`$ $`=`$ $`\sqrt{2\pi }\sigma ^2q(E)\mathrm{e}^{\sigma ^2q^2(E)/2}.`$ (59)
The zeroth order approximation, $`g_E(E)=\delta (Eϵ_{max})`$, Eq.(34), is largely modified at one loop. The imaginary part of the self-energy vanishes at $`0`$ and $`ϵ_{max}`$ since both $`𝒢(q(0),p)`$ and $`𝒢(q(ϵ_{max}),p)`$ are zero. The zeroth order delta function is moved to a value $`E^{}`$ which verifies the relation :

$$E^{}=ϵ_{max}+Re\mathrm{\Sigma }_1(\mathrm{},E^{}),$$
(60)
and the density of states at this order is

$`g_E(E)`$ $`=`$ $`\sqrt{{\displaystyle \frac{2\sigma ^2}{\pi }}}{\displaystyle \frac{q(E)\mathrm{e}^{\sigma ^2q(E)^2/2}}{\left(E\mathrm{\Sigma }_1(\mathrm{},E)\right)^2+q(E)^2\mathrm{e}^{\sigma ^2q(E)^2}}},0Eϵ_{max},`$ (61)
$`g_E(E)`$ $`=`$ $`{\displaystyle \frac{1}{\left|\frac{d\left[Re\mathrm{\Sigma }_1(\mathrm{},E)\right]}{dE}|_{E=E^{}}\right|}}\delta (EE^{}),E>ϵ_{\mathrm{max}}.`$ (62)
### C Two loops

The order 1/$`\rho ^2`$ arises from including $`_{R2}^{(R)}(p)`$ on the geometric series for the resolvent. The particle label repetitions can basically arise in four ways, that we schematically depict as

$`\mathrm{}1\mathrm{}2\mathrm{}2\mathrm{}1\mathrm{}`$ (63)
$`\mathrm{}1\mathrm{}2\mathrm{}1\mathrm{}2\mathrm{}`$ (64)
$`\mathrm{}1\mathrm{}1\mathrm{}1\mathrm{}`$ (65)
$`\mathrm{}1\mathrm{}1\mathrm{}2\mathrm{}2\mathrm{}`$ (66)
The repetitions depicted in (63) and (64) are genuine contributions to the self-energy. The repetitions (66) are contributions to the Dyson resumation of $`\mathrm{\Sigma }_1`$, (the $`G_0^3(p,z)\mathrm{\Sigma }_1^2(p,z)`$ term in Eq.(39)), while the repetitions (65) contribute both to $`\mathrm{\Sigma }_2`$ and to $`G_0^3(p,z)\mathrm{\Sigma }_1^2(p,z)`$.

In Appendix A we list all the two-loops terms giving the contribution to the self energy $`\mathrm{\Sigma }_2`$. Because of the complexity of the result, we shall focus herafter only onto the low momentum regime. For $`p=0`$, there is a cancellation among all the various diagrams and the self energy vanishes as it should. The leading contribution to the imaginary part of the self-energy scales as

$`Im\mathrm{\Sigma }_2(p,E)=Cp^2q^{d2}(E)`$ (67)
In the low energy, low momentum region, close to the peak $`pE^{1/2}`$ of the structure factor, it becomes:

$$\mathrm{\Sigma }_2p^2E^{(d2)/2}p^d.$$
(68)
Consequently, for small values of the external momentum $`p`$, the contribution to the width arising from the $`1/\rho ^2`$ order is much broader than the contribution from the $`1/\rho `$ order. In the frequency domain, the expression (68) gives indeed the broadening:

$$\mathrm{\Gamma }p^{d1}$$
(69)
i.e. $`p^2`$ in the three dimensional case, very different from the result (48).

Putting all together, our final expression for the imaginary part of the self-energy in the low momentum regime, for a function $`f(r)=\widehat{f}(r/\sigma )`$, is up to second order:

$$Im\mathrm{\Sigma }=\rho \stackrel{~}{f}(0)\left[\frac{1}{\rho \sigma ^d}(p\sigma )^{d+2}A_1+\left(\frac{1}{\rho \sigma ^d}\right)^2(p\sigma )^dA_2\right]$$
(70)
where the coefficients $`A_1`$ and $`A_2`$ are pure numbers $`𝒪(1)`$ which can be determined after choosing the adimensional function $`f(x)`$ (the gaussian case will be studied in the next section). From the qualitative point of view it is worthwhile to note that (70) describes the most general case. It is quite easy to realize that even considering other multi-loops terms, for a small momentum, the leading order term would be of the form (69).

Let us note that the parameter of such expansions is actually the inverse of the number of particles $`N_R\rho \sigma ^d`$, lying inside a cube whose dimensions are given by the range of the interaction. Let us picture the whole scenario as it emerges from two-loops computation:

* As a very general feature of topologically disordered systems, we obtain the broadening $`p^{d1}`$ of the resonance peak of the $`S(p,\omega )`$ at low enough momentum $`p`$, where the two-loop contributions is dominant. Because of the $`p^{d1}`$ behaviour, topologically disordered systems spread the energy of an incoming plane wave on a wider range of frequencies than non-topologically disordered systems. This is a quantitative characterization of the intuitive impression that they are “more disordered”. However, this feature is rather subtle, and only appears when the contributions $`1/N_R^2`$ to the self energy are considered.
* At very high densities there exists a characteristic values $`p_c`$ where there is a cross-over from the $`p^{d1}`$ regime to the $`p^{d+1}`$ regime. From (70) the cross-over momentum is easily found:

$$(\sigma p_c)^2=\frac{A2}{A1}\frac{1}{\rho \sigma ^3}\frac{1}{N_R}$$
(71)
* For $`p>p_c`$ the one loop term mainly contributes and the broadening is controlled by $`\mathrm{\Sigma }_1`$. In this regime, the asymptotic result (48) has to be taken with a lot of care. In fact it is supposed to be correct only for $`\sigma p<<1`$, a regime which basically is dominated by the two-loops contributions. If $`\sigma p_c`$ is not small enough to be still in the same regime we are no longer legitimate to assume that (48) give the correct results and we have to consider the full one-loop result (40).

## III Comparison with numerical simulations

The dynamical structure factor $`S_E(q,E)`$ can be numerically computed. Its definition, Eq.(6), implies that it is a well normalized distribution function,

$$S_E(p,E)0,_{\mathrm{}}^{\mathrm{}}dES_E(p,E)=1.$$
(72)
The moments of this distribution function can be calculated in terms of the powers of the matrix $`M`$ defined in Eq.(5),

$$_{\mathrm{}}^{\mathrm{}}𝑑ES_E(p,E)E^R=\overline{\frac{1}{N}\underset{i,j}{}\mathrm{e}^{\mathrm{i}p(x_ix_j)}\left(M^R\right)_{i,j}}.$$
(73)
If the function $`f`$ is non-negative, one can cut the previous integrals at $`E=0`$, since the quadratic form associated to the matrix $`M`$ is semi positive-definite:

$$\underset{i,j}{}\phi _iM_{i,j}\phi _i=\frac{1}{2}\underset{i,j}{}f(x_ix_j)(\phi _i\phi _j)^20$$
(74)
Similarly, one can obtain the density of states replacing Eq.(73) by

$$_{\mathrm{}}^{\mathrm{}}𝑑Eg_E(E)E^R=\overline{\frac{1}{N}\underset{i,j}{}\frac{v_iv_j}{_kv_k^2}\left(M^R\right)_{i,j}},$$
(75)
where the $`v_i`$ are random numbers chosen with uniform probability between $`1`$ and $`1`$, and the overline now means average over the particle positions and the $`v_i`$. Thus we see that one can uses the method of moments  complemented with a truncation procedure . Using the gaussian function $`f`$ of Eq.(49) truncated at four $`\sigma `$’s, we have been able to reconstruct the $`S_E(p,E)`$ on systems with up to 32768 ($`32^3`$) particles, using 100 moments. All the simulations in this paper have been done generating ten samples of the disordered configurations: on a box of side $`L=32\sigma `$, we place at random $`\rho L^3`$ particles. We apply periodic boundary conditions, the minimum available momentum thus being $`\frac{2\pi }{32\sigma }`$. Regarding the choice of the density, let us recall that our natural length-unit is $`\sigma `$, that can be considered as the analogous of the wavelength corresponding to the first maximum of the dynamical structure factor. It is therefore clear that the densities that correspond with the experimental situation are of the order of $`\sigma ^3`$. In this paper we have explored the range $`0.2\sigma ^3\rho 1\sigma ^3`$. In what follows, momentum will be measured in units of $`\sigma ^1`$ and density in units of $`\sigma ^3`$.

The choice of the number of calculated moments is conditioned by two conflicting goals. On the first place, the irregularities on the reconstructed probability distribution grows significantly with the number of calculated momenta, which is due either to the statistical fluctuations on the finite number of generated samples, or to round-off errors on the Lanczos recursion. On the other hand using a too small number of moments can smooth-out real features of the curves. As a rule we show the curve reconstructed with the minimum number of moments for which the width of the $`S_E(p,E)`$ is stable. We should emphasize that the moments method is a statistical one. With a finite number of particles, there are only $`\rho L^3`$ available eigenvalues. If the number of eigenvalues around the maximum of the $`S_E(p,E)`$ is small, the method of moments does not yield meaningful results. Moreover, since the $`S_E(p,E)`$ is an statistical quantity itself, meaningful results cannot be obtained even with a full diagonalization of the matrix $`M`$. In order to get a feeling on what difficulties might be encountered, it is useful to calculate the first and second moments of the $`S(p,E)`$:

$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑ES_E(p,E)E`$ $`=`$ $`ϵ(p),`$ (76)
$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑ES_E(p,E)E^2`$ $`=`$ $`ϵ(p)^2+{\displaystyle \frac{2\rho }{(2\pi )^d}}{\displaystyle 𝑑k\stackrel{~}{f}(k)[\stackrel{~}{f}(k)\stackrel{~}{f}(pk)]}.`$ (77)
It is clear in what sense in the $`\rho \mathrm{}`$ limit the $`S_E(p,E)`$ becomes a Dirac delta: its mean-value, $`ϵ(p)`$, grows like $`\rho `$, while its variance only grows like $`\sqrt{\rho }`$. In particular, one can take the $`p\mathrm{}`$ limit in the above equations, and see what happens with the density of available eigenvalues. For the model described in Eq. (49) , one obtains:

$`{\displaystyle _0^{\mathrm{}}}𝑑Eg_E(E)E`$ $`=`$ $`\rho (2\pi \sigma ^2)^{3/2},`$ (78)
$`{\displaystyle _0^{\mathrm{}}}𝑑Eg_E(E)E^2`$ $`=`$ $`\rho ^2(2\pi \sigma ^2)^3+2\rho (\pi \sigma ^2)^{3/2}.`$ (79)
The above results are rather discouraging, because one immediately observes that at fixed $`p`$, the larger is $`\rho `$, the lesser states are around $`ϵ(p)`$. Figure 1 may help in clarifying this point.

The statistical sampling improves significantly with decreasing densities, where however our analytical calculation is less likely to work. In practice, the above considerations imply that the smaller is $`\rho `$, the smaller is the minimum momentum for which the $`S_E(p,E)`$ can be safely estimated on a finite box.

As a first comparison between the analytical and the numerical computations, in figure (2) we show the DOS obtained numerically and the (continuous part) of the DOS at one-loop level as obtained by means of the eq. (62) (neglecting the contribution of the delta function). As it has been already pointed out in , the high density expansion is not really suitable to describe the DOS in the neighborhood of the maximum. In fact, we have shown that at one loop the zeroth order delta function splits from the continuous part of the DOS, carrying still a finite weight. Furthermore the continuous part (62) is defined only for $`E<ϵ_{max}`$; in order to correctly compute the density of states in this region we need at least a partial resummation of the higher orders in the $`\rho ^1`$ expansion, in the same way that it was done in  for non translational invariant potentials. This task is not out of reach, however it goes beyond the scope of this paper. Because of this limitations, here only the low eigenvalues region of the spectrum can be safely compared. We see that the numerical curve obtained with a large number of moments is badly oscillating, because the single delta functions which form the spectrum of a finite matrix are correctly reproduced by the method of moments. If we consider the curve obtained considering only $`10`$ moments instead, there is a smoothing of the previous curves and we see that it is in very nice agreement with the one loop computation.

Let us turn to the comparison of the dynamical structure factor. One would like to check wether all the scaling regimes identified in subsection II C really hold. One should therefore know at which value of the exchanged momentum the one loop result will start to be dominating. The two coefficients of Eq.(70), for the Gaussian choice (49), in the $`3d`$ case the coefficients turn out to be:

$`A_1`$ $`=`$ $`7/12\mathrm{\Omega }_3`$ (80)
$`A_2`$ $``$ $`1.25\mathrm{\Omega }_3^2`$ (81)
with

$$\mathrm{\Omega }_3=\frac{d\mathrm{\Omega }_3}{(2\pi )^3}$$
(82)
The crossover momentum can then be estimated to be $`0.31/\sigma `$ for $`\rho =1\sigma ^3`$ and $`0.71/\sigma `$ for $`\rho =0.2\sigma ^3`$. In a box of size $`L=32\sigma `$ the first available momentum is $`0.19/\sigma `$, but unfortunately this is not the first momentum for which one can reconstruct accurately the $`S_E(p,E)`$. Therefore, we should expect that our data can be accurately described using only by the full one-loop expressions. Of course, one could try to compare with still lower densities, but this would be both risky (ours is a $`1/\rho `$ expansion) and unphysical (glasses have densities of order one, if measured in the natural length units of their interaction potentials).

The first comparison is for the density $`\rho =1`$. As can be seen from Fig. 2, for a box with $`32^3`$ particles, when the position of the peak is at $`E_{peak}<12`$ there are too few states and the numerical results have no statistical significance. Hence, the lowest value of the external momentum which is possible to compare safely with the theory turns out to be $`p=0.95/\sigma `$ (i.e. $`n=5`$). The corresponding results are displayed in Fig. 3. We can see that for $`n=5,7`$ the one-loop computation is in very good agreement with the numerical data. When the momentum is too large ($`n=9,11`$ in the figure), $`E_{peak}`$ becomes comparable to $`ϵ_{max}`$ and the agreement cannot be longer so good. The analytical curve indeed is defined only up to $`ϵ_{max}`$ (at all orders in perturbation theory) and the tail which follows the peak cannot be reproduced by the theory without further developments.

The next comparison is at density $`\rho =0.2`$. Even if one might expect that this density is too low for the $`1/\rho `$ expansion to be accurate, the one-loop computation still describes quite well the numerical curves, up to the region where $`E_{peak}`$ and $`ϵ_{max}`$ become comparable. One is tempted to conclude that the difference in fig. (4) between theory and numerics for the curves which are not strongly affected by the tail effect ($`n=4,6`$) are due mainly to the following terms of the series in the $`1/\rho `$ expansion. In particular the correction of the broadening is of order $`10\%`$ while the position of the peak seems to be less affected by higher order corrections, specially for the $`n=4`$ curve.

As we have pointed out in the introduction, a controversial point is the existence of a scaling law for the broadening of the spectral peak with respect to $`p`$, in the range of momenta explored by X-rays scattering experiments. On the other hand, there is a general agreement on the existence of the dispersion law for the peak position of the kind $`E_{peak}p^2`$ ($`\omega ^{peak}p`$), up to momenta where the corresponding wavelength is comparable with the mean separation of the particles, which enable us to identify the peak as an acoustic one. Having shown that our numerical data are reasonably described by the one loop computation, let us see what are the predictions for these spectral features at this level of accuracy. In the figure (5) we plot the real and the imaginary part of $`\mathrm{\Sigma }_1`$ at the point $`E_{peak}`$, obtained by means of (56,57) as functions of $`p`$, for the density $`\rho =1`$

The real part follows the dispersion law $`p^2`$ up to momentum of order $`1/\sigma `$, and no significant dependence on the density is observed while that scaling law still holds, see Eqs.(46) and (47). For the imaginary part, we already know from the asymptotic analysis of the $`p0`$ regime, that the one-loop low momentum behaviour is $`p^5`$ ($`p^4`$ in the frequency domain), however the range of momenta which are studied by numerical simulations, (and by IXS scattering experiments!) surely are not in that asymptotic regime. The figure shows clearly that the lowest momentum that we are able to numerically study is already in a region where that scaling law does not hold anymore. As a matter of fact, the embarrassing conclusion is reached that the scaling of the width of the dynamical structure factor at one loop (which coincides with the numerical simulation) could mimick a $`p^2`$ behaviour if the range of momenta which is studied is not wide enough.

## IV About Particle Correlations

The formula (17) is a good starting point to speculate about the possibility of extending our computation to the case where the particles are not chosen with a random distribution but, for example, with the equilibrium Gibbs distribution. In that case, the average on the position of the particles is made using :

$$\frac{1}{V^{R+1}}\underset{i}{\overset{R+1}{}}d^dx_ig^{(R+1)}(x_1\mathrm{}x_{R+1}),$$
(83)
where $`g^{(R+1)}`$ is the $`R+1`$-points correlation function. Although the computation using the full correlation function would be exceedingly difficult, some progress can be made by using the so called superposition approximation:

$$g(x_1\mathrm{}x_{R+1})=g(x_1x_2)g(x_2x_3)\mathrm{}g(x_Rx_{R+1})$$
(84)
where the pair correlation function is used to take into account the correlation of the position of the particles. The superposition approximation has been probed to be reasonably correct in describing the high order static correlation functions of the supercooled liquids, both in the computation of the coupling coefficients in the mode coupling theory  and in the computation of the vibrational entropy .

It is not difficult to see that, to the lowest order in $`1/\rho `$, the superposition approximation can be embebed in our calculation if we make the substitution:

$$f(r)g(r)f(r)$$
(85)
This is rather important, because, for typical applications, the function $`f`$ being badly divergent at low distances, does not a have a Fourier transform. On the other hand, the function $`g(r)`$ typically tends to zero at the origin exponentially, thus taking care of the algebraic divergence of $`f(r)`$.

However, the substitution (85) cannot be blindly applied to the calculation of the $`1/\rho `$ corrections. Let us consider for instance the second moment (77), that in the correlated case yield (let $`F`$ be the Fourier transform of f(r)g(r))

$`\widehat{ϵ}(p)`$ $`=`$ $`\rho F(0)\rho F(p),`$ (86)
$`\overline{M^2(p)}`$ $`=`$ $`[\widehat{ϵ}(p)]^2+2\rho {\displaystyle 𝑑xg(x)f^2(x)\left(1e^{ipx}\right)}`$ (87)
while the subsitution (85) would have produced a factor $`g^2(x)f^2(x)`$ in the above integral. A little thought reveals that, for the $`1/\rho `$ calculation reported in Eq.(36), this problem will arise whenever there is a string of only diagonal terms ($`\delta _{k_i,k_{i+1}}f(x_{k_i}x_{z_i})`$) between the two repeated particle labels. The diagonal terms are simply $`[\rho F(0)]^b/z^{b+1}`$ that add up to

$$\frac{1}{z\rho F(0)}.$$

Therefore, we conclude that the self-energy at order $`1/\rho `$ in the superposition approximation is

$`\mathrm{\Sigma }_1(p,z)`$ $`=`$ $`{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{dq}{(2\pi )^D}[\rho F(pq)\rho F(q)]^2\left(\frac{1}{z\widehat{ϵ}(q)}\frac{1}{z\rho F(0)}\right)}+`$ (88)
$`+`$ $`2\rho {\displaystyle 𝑑xg(x)f^2(x)\frac{\left(1e^{ipx}\right)}{z\rho F(0)}}.`$ (89)
It is clear that none of the two new terms contribute to the imaginary part close to $`\widehat{ϵ}(p)`$, and therefore our discussion of the broadening of the spectral line, at this order of the $`1/\rho `$ expansion is basically unchanged. One could similarly discuss the peculiarities of the $`1/\rho ^2`$ term, but this is left for future work. However, the single term

$$2\rho 𝑑xg(x)f^2(x)\frac{\left(1e^{ipx}\right)}{z\rho F(0)},$$
(90)
points to the possibility of a qualitatively new behaviour. Indeed, this term is the first of a full category, in which the repeated particle labels always corresponds to the same pair of particles, separated by a string of diagonal terms. It is not hard to show that the contribution to the self-energy of this family of terms is:

$$\frac{\rho }{2}\underset{n=1}{\overset{\mathrm{}}{}}𝑑xg(x)\left(1e^{ipx}\right)\frac{[2f(x)]^{n+1}}{[z\rho F(0)]^n}=2\rho 𝑑xg(x)f^2(x)\frac{\left(1e^{ipx}\right)}{z\rho F(0)2f(x)}.$$
(91)
This result implies that when the function $`f(x)`$ diverge at small distance, the density of state will have a tail to infinity, exponentially supressed by the $`g(x)`$ function. This exponential tail is of a different origin from the instantonic contribution recently calculated in Ref.. Indeed the instantonic exponential tail arise also for a fully regular f(x), due to anomalously dense regions.

## V Discussion, Conclusions and Outlooks

The aim of this work was to study the spectral properties of euclidean random matrices using an expansion on the inverse particle density number. That approach is supposed to describe correctly topologically disordered systems, in the context of the Instantaneous Normal Modes approach to their high-frequency dynamics. In particular, we have performed a two-loop computation of the dynamical structure factor, where the expansion parameter is the inverse of the density of particles. The results have been analytically studied in the low momentum regime, where the disorder reveals itself by the broadening of the resonance peak. The width of that spectral line has been computed up to order $`1/N_R^2`$, where $`N_R`$ is the number of particles which effectively interact with a given particle. The behaviour of the width $`\mathrm{\Gamma }`$ in the frequency domain, as function of the momentum $`p`$ has been shown to follow a quite complicated law. When the exchanged momentum, $`p`$, is quite smaller that the inverse wavelength of the first maximum of the static structure factor, the width of the spectral line follows $`\mathrm{\Gamma }p^2`$. This seem to be an intrinsic difference between topologically-ordered and topologically-disordered random systems, the Rayleigh $`p^4`$ scaling being always found for the formers. At a momentum depending on the particle density and the details of the potential, the $`1/N_R`$ contribution to the self-energy becomes dominant, and the $`p^4`$ scaling appears. Finally, when the exchanged momentum is of the order of the inverse wavelength of the first maximum of the static structure factor, the dynamical structure factor starts its collapse onto the density of states. The width of the spectral line cannot grow at such a violent pace as $`p^4`$, and if a not too large momentum-range is explored, it can certainly mimick a $`p^2`$ scaling. It is possible that the $`p^2`$ scaling found in the IXS experiments correspond to the later transient behaviour. Nevertheless, it would be very interesting to explore experimentally a wider range of exchanged momentum (maybe combining several scattering probes), in order to test this scaling picture with three separated regimes.

It would be also very useful to compare the predictions of this computation with the ones which have been obtained in the framework of mode coupling approximation . There an hard sphere system is investigated and a $`p^2`$ law is obtained only for very low values of the frequency of the peak, which at the moment are not accessible to the experiments, while a crossover to a quite different behaviour is present when the external momentum moves to higher values, comparable with the ones explored by the IXS scattering experiments. Even in that computation no clear scaling law actually emerges and qualitatively that scenario is quite similar to the one emerging from our approach. A quantitative comparison though, is not straightforward because our analytic computation does not concern any particular model of topologically disordered system, the choice (49) being quite generic in order to understand which are the main features which determine the particular kind of broadening of the resonance peak.

The main point, however, is the reliability of our computation in understanding the experimental findings and, possibly, in predicting features not yet observed. Apparently, if one wants to apply the computation presented here to the study of the spectral properties of realistic systems, such as silica or fragile glasses, one faces the following limitations:

* The computation has been explicitly performed in the case where particles are uncorrelated which is deeply different from the situation found in glasses and supercooled liquids.
* The values of the coefficients $`A_1`$ and $`A_2`$ in (70) seems to depend very badly on the choice of the interaction $`f(r)`$ hence at this level is quite unclear the generality of the predictions that one could obtain. from (71) before performing the appropriate computations for a given system.

However, in section IV it has been shown that the correlations between the particles can be at least partially taken into account by performing the transformation (85), which amounts to dress the bare interaction $`f(r)`$ between two particles at a given distance $`r`$, with the probability, given by the static structure factor $`g(r)`$ of the real system to be studied, to find them at that distance $`r`$. Furthermore this procedure corresponds to the well-known superposition approximation of the correlation functions of the theory of liquids, a kind of approximation that has been largely exploited previously in different analytical approaches to the study of the glassy phase. More interestingly, the superposition approximation makes even the latter limitation less severe. As a matter of fact the Gaussian choice (49) has to mimic not the bare interaction $`f(r)`$, that would be an hopeless task, but the function $`f(r)g(r)`$ which sounds more feasible. We have also shown how the terms for which the substitution $`f(r)f(r)g(r)`$ cannot be made, give rise to an exponential tail of the density of states, of non instantonic (see ) origin. More in detail, by considering only the bare interaction one should take:

$$f(r)=\frac{d^2V(r)}{dr^2}$$
(92)
where $`V(r)`$ is the pair potential of the particles, which typically has a strong repulsive barrier at $`r0`$, missing in the Gaussian choice. The static structure factor $`g(r)`$, on the other hand, is $`e^{\beta V(r)}`$ at $`r0`$, hence the function $`f(r)g(r)`$ is typically formed by a single peak in the position $`r_0`$, where the static structure factor has its maximum, is not singular at $`r0`$ and it does have a finite Fourier transform looking quite similar to (50). As a first approximation, the above argument yields $`\sigma r_0`$, i.e. a ’few’ Angstrom for realistic systems, implying that X-rays experiments surely cannot explore the low momentum region but they can observe the region $`p\sigma >0.20.3`$, which is roughly the same region as that investigated in the numerical simulations.

Hence we expect that the computation we have presented here, if suitably modified for including the details of the system under study, will allow to compute the dynamical structure factor of structural glasses and supercooled liquids.

## Acknowledgments

We would like to thanks G.Biroli, W.Kob, R.Monasson, G.Ruocco and G.Viliani for valuable discussions and crucial suggestions. This work has been performed in part during the stay of two of us (MM and PV) at the Institute for Theoretical Physics of the University of California, Santa Barbara. We thank the Institute for its hospitality and acknowledge the partial support by National Science Foundation under grant No.PHY94-07194. V.M.-M. is a M.E.C. fellow. P.V. would like also to thank the SPHYNX program, which has supported his staying at the Ecole Normale Superieure in Paris.

## APPENDIX A

In this appendix we shall show the two loops integrals and the corresponding diagrams arising from the $`1/\rho `$ expansion of the resolvent to order 1/$`\rho ^2`$. It turns out to be useful to express the integrals in terms of the propagators $`\rho \stackrel{~}{f}(q)`$ and $`G_s(q)`$, where

$`G_s(q)a^2\left(G_0(q){\displaystyle \frac{1}{a}}\right)`$ (93)


text:
# Response to Fackerell’s Article

## 1 Introduction

In an article titled “Remarks on the Yılmaz and Alley Papers” E.D. Fackerell  claims that the Yılmaz metric, used by H. Yılmaz  and C.O. Alley  for the treatment of parallel slabs in general relativity, is not a solution to the field equations of general relativity. His train of thought in making this statement seems to be as follows:

The Yılmaz metric used in the said treatments is of the form

$$ds^2=e^{2\varphi }dt^2e^{2\varphi }[e^{2ϵ\varphi }(dx^2+dy^2)+e^{4ϵ\varphi }dz^2]$$
(1)
$$\varphi =az+\frac{1}{2}\sigma z^2+C$$
(2)
where $`ϵ=\pm 1`$ for general relativity. \[$`ϵ=0`$ for the Yılmaz theory.\] But Fackerell, based on some calculations of his, seems to conclude that the field equations will everywhere require $`g_{zz}`$ to be

$$g_{zz}=e^{2\left(1+2ϵ\right)\varphi +2\mathrm{ln}\left|_z\varphi \right|+K}.$$
(3)
Therefore, he says, the Yılmaz metric will be valid if and only if $`_z\varphi =`$ Constant so that one can set $`2\mathrm{ln}|_z\varphi |+K=0`$ to regain the Yılmaz form. Since $`_z\varphi `$ is not a constant in the Yılmaz field $`\varphi `$ in Eq. (2) above, the metric used in the treatment, he claims, cannot be a solution to the field equations of general relativity!

However, a simple calculation, which can be done by hand, shows that the Yılmaz metrics for $`ϵ=\pm 1`$ satisfy the field equations of general relativity everywhere exactly. If Fackerell did this calculation he would have seen that it is indeed true. Therefore the question is not whether Yılmaz and Alley used wrong metrics in their works but rather how and where precisely Fackerell himself made a fundamental error in his calculations.

## 2 Fackerell’s Fundamental Error

To find this out, we shall first verify that Yılmaz’ metrics are indeed solutions of the required field equations. Fackerell’s error will be clearly seen in the course of this exercise. The general relativity equations to satisfy with $`\tau _\mu ^\nu `$ as a matter tensor are

$$\frac{1}{2}G_\mu ^\nu =\tau _\mu ^\nu .$$
(4)
There are two distinct cases with two sub-cases each: \[Note 1\]

$$\mathrm{A}):\sigma =0,ϵ=\pm 1(\mathrm{Matter}\mathrm{free}\mathrm{regions})$$

$$\mathrm{B}):\sigma 0,ϵ=\pm 1(\mathrm{Regions}\mathrm{with}\mathrm{matter})$$

where $`\sigma `$ corresponds to the ordinary Laplacian of $`\varphi `$. The Yılmaz metric is of the form

$$ds^2=e^{2\varphi }dt^2e^\mu (dx^2+dy^2)e^\eta dz^2.$$
(5)
where $`\varphi `$, $`\mu `$ and $`\eta `$ are functions only of the coordinate z. Using the shorthand $`(^{})=_z\varphi `$, the general relativity equations (4) for metric (5), after multiplying by $`8e^\eta `$, are:

$$\mu ^{}(3\mu ^{}2\eta ^{})+4\mu ^{\prime \prime }=8e^\eta \tau _0^0$$
(6)
$$4\varphi ^{\mathrm{\hspace{0.17em}2}}+(\mu ^{}2\varphi ^{})(\mu ^{}\eta ^{})+2(\mu ^{\prime \prime }2\varphi ^{\prime \prime })=8e^\eta \tau _1^1=8e^\eta \tau _2^2$$
(7)
$$\mu ^{}(\mu ^{}4\varphi ^{})=8e^\eta \tau _3^3.$$
(8)
Let us first consider the matter-free case (A):
For the Yilmaz metric, where $`\mu =2(1+ϵ)\varphi `$ and $`\eta =2(1+2ϵ)\varphi `$, all terms quadratic in $`\varphi ^{}`$ cancel. (The algebra is exhibited in \[Note 2\].) The equations (6) and (7) in this case become:

$$(1+ϵ)\varphi ^{\prime \prime }=0$$
(9)
$$\frac{1}{2}ϵ\varphi ^{\prime \prime }=0$$
(10)
with the solution

$$\varphi ^{\prime \prime }=0$$
(11)
$$\varphi =az+C$$
(12)
where $`a`$ and $`C`$ are constants, hence $`\varphi ^{}=`$Constant is automatic. Setting $`2\mathrm{ln}|a|+K=0`$ in equation (3) above, as Fackerell himself suggests, we have

$$g_{zz}=e^{2\left(1+2ϵ\right)\varphi }$$
(13)
as originally stated by Yılmaz.
Let us now consider also the regions with matter, case (B):
The $`\tau _\mu ^\nu 0`$ equations (6) - (8) are satisfied with the matter tensor of the form

$$\tau _\mu ^\nu =e^\eta \mathrm{diag}[(1+ϵ,ϵ/2,ϵ/2,0)\sigma ]$$
(14)
where by directly substituting $`\mu =2(1+ϵ)\varphi `$, $`\eta =2(1+2ϵ)\varphi `$, all terms quadratic in $`\varphi ^{}`$ again cancel for $`ϵ=\pm 1`$ (as shown in \[Note 2\]) and we have

$$(1+ϵ)\varphi ^{\prime \prime }=(1+ϵ)\sigma $$
(15)
$$\frac{1}{2}ϵ\varphi ^{\prime \prime }=\frac{1}{2}ϵ\sigma $$
(16)
with the solution

$$\varphi ^{\prime \prime }=\sigma $$
(17)
$$\varphi =az+\frac{1}{2}\sigma z^2+C$$
(18)
as stated by Alley and Yılmaz. But now there are no logarithmic terms in the exponential because $`\varphi ^{\prime \prime }`$ terms causing them are taken away by equating them to the matter terms. The result is again the form of equation (13) above

$$g_{zz}=e^{2\left(1+2ϵ\right)\varphi }$$

as originally stated by Yılmaz.

The whole discussion can be reduced into the realization that

$$\varphi ^{\prime \prime }=\{\begin{array}{cc}0& (\mathrm{Outside}\mathrm{a}\mathrm{slab}),\\ \sigma & (\mathrm{Inside}\mathrm{a}\mathrm{slab})\end{array}$$
(19)
Nothing else is needed, not even Fackerell’s way of removing the logarithmic terms above, since the form given by Yılmaz makes it superfluous. \[Note 3\]

Thus, whether we like it or not, the Yılmaz metric is a solution of the field equations of general relativity, both in the matter-free regions and in the regions with matter. Then, where did Fackerell make his mistake? In the light of the above discussion we do not have far to look. Since in a matter-free region $`\varphi `$ is linear $`(\varphi ^{\prime \prime }=0`$) Fackerell thought this would contradict the quadratic term in equation (18) above

$$\varphi =az+\frac{1}{2}\sigma z^2+C.$$

But those are different regions with different material contents as indicated above. Linearity in coordinates of a potential in one region does not prevent it from being nonlinear in coordinates in a different region. Thus, contrary to this fact, Fackerell is requiring that the linearity of $`\varphi `$ in the matter-free region be imposed on the region with matter. This is equivalent to requiring that there be no matter terms at all, anywhere. A more elementary calculational error is hardly imaginable.

We should really stop here and go no further with Fackerell’s paper. All his other asides such as his appeals to authority and his prima facie case, etc., are irrelevant since his main tenet is false. Therefore it is useless to go over every conceivable argument. We must however, mention just one additional error, because this one is so hopelessly hidden that the reader, not being able to penetrate, might think he has something.

## 3 Fackerell’s Second Error: His Treatment of the Parallel Slab Problem in General Relativity

Let us now also examine Fackerell’s own treatment of the parallel slabs. Fackerell claims that the appropriate metric to use for the slab calculations is the Taub metric

$$ds^2=f^{1/2}dt^2f(dx^2+dy^2)f^{1/2}dz^2$$
(20)
where

$$f=(1+kz).$$
(21)
However, he does not use the slabs as originally defined, and he does not use equations of motion. Instead, he considers infinitely thin sheets. Under some sophisticated looking maneuvers he claims he gets certain results which sound impressive. But it can be shown by a simple calculation that, for the Taub metric the left-hand side (hence also the right-hand side) of the field equations are identically zero. Again, if Fackerell performed this calculation he would have seen that it is true. In any case the proof is so simple that we shall here reproduce it: First, identify from the Taub metric by comparing with the general metric, equation (5), that $`\mu =4\varphi `$, $`\eta =2\varphi `$. Then simply calculate to get,

$$e^{4\varphi }=f=1+kz$$
(22)
$$\varphi =\frac{1}{4}\mathrm{ln}(1+kz)$$
(23)
$$\varphi ^{}=(k/4)/(1+kz)$$
(24)
$$\varphi ^{\prime \prime }=(k^2/4)/(1+kz)^2=4\varphi ^{\mathrm{\hspace{0.17em}2}}.$$
(25)
Finally, substituting into the required field equations (6) - (8) one finds that all expressions on the left-hand side are identically zero. (The steps in this calculation are exhibited in \[Note 4\].) Therefore the right-hand side is

$$\tau _\mu ^\nu =\mathrm{diag}[(0,0,0,0)].$$
(26)
This is because, unlike the Yilmaz metrics, the term linear in $`\varphi ^{\prime \prime }`$ (that is, the Laplacian) also disappears and nothing is left to represent matter.

As to the equations of motion, we shall then have

$$\sigma d^2x_\mu /ds^2=\frac{1}{2}_\mu g_{\alpha \beta }\tau ^{\alpha \beta }=0.$$
(27)
This is clearly the answer to the problem of motion with the Taub metric. But as we have just said, Fackerell does not use the equations of motion of the theory, and he does not use the definition of the slab as originally intended. Instead, he introduces infinitely thin sheets. This is presumably to get around the difficulty that he cannot introduce a matter density $`\tau _\mu ^\nu `$. The truth is that as long as he uses the Taub metric he cannot introduce any matter at all. Under these circumstances how can the “junction conditions” etc., be of any help? Can there be a magic in the junction conditions so as to create something out of absolutely nothing? \[Note 5\]

Yılmaz and Alley do not use such devices. When they say slab, they mean slab, when they say equations of motion, they mean equations of motion. They do not use junction conditions to supplant the equations of motion nor do they stray from the boundary conditions of the standard field theory. As Fackerell himself points out, the “junction conditions” he uses are not necessarily equivalent to the boundary conditions of the standard field theory. Fackerell’s deductions from the Taub metric are clearly false.

## 4 Discussion

We have been informed that during a recent discussion on the internet newsgroup $`\mathrm{𝑠𝑐𝑖}.\mathrm{𝑝ℎ𝑦𝑠𝑖𝑐𝑠}.\mathrm{𝑟𝑒𝑙𝑎𝑡𝑖𝑣𝑖𝑡𝑦}`$, a respected member of the group stated:
“Alley and Yılmaz’s claims along these lines have been demolished, thoroughly and in detail, in a paper by Edward Fackerell, …\[ref.  of this paper\]…, Fackerell shows that Yılmaz and Alley made a careless mistake in their calculation, and he gives the correct calculation, for both flat plates and spherical shells.”

We cannot expect everyone to dig into the truth of the matter but we must urge Fackerell to honorably retract his position. This sort of thing is not doing any service to science – nor any justice to the scientific truth.

Acknowledgment
We wish to thank Kirk Burrows for carrying out the detailed calculations associated with our treatment of the two parallel slabs problem . His expert knowledge of the symbolic calculational capabilities of Mathematica has been invaluable in our continuing investigations of the theories of Yılmaz and of Einstein.

Notes

\[Note 1\] Fackerell treats the cases $`ϵ=+1`$ and $`ϵ=1`$ separately from the start. This gets confusing and repetitive. We treat them uniformly with $`ϵ`$ as a parameter and only at the end set $`ϵ=+1`$ and $`ϵ=1`$.

\[Note 2\] Substituting the expressions $`\mu =2(1+ϵ)\varphi `$ and $`\eta =2(1+2ϵ)\varphi `$ into the left-hand sides of equations (6) - (8), one obtains the following:

Left-hand side of Equation (6):

$$\mu ^{}(3\mu ^{}2\eta ^{})+4\mu ^{\prime \prime },$$

$$2(1+ϵ)\varphi ^{}[6(1+ϵ)\varphi ^{}4(1+2ϵ)\varphi ^{}]+4[2(1+ϵ)\varphi ^{\prime \prime }],$$

$$2(1+ϵ)\varphi ^{}[2(1ϵ)\varphi ^{}]+8(1+ϵ)\varphi ^{\prime \prime },$$

$$4\varphi ^{\mathrm{\hspace{0.17em}2}}(1ϵ^2)+8(1+ϵ)\varphi ^{\prime \prime }$$

Left-hand side of Equation (7):

$$4\varphi ^{\mathrm{\hspace{0.17em}2}}+(\mu ^{}2\varphi ^{})(\mu ^{}\eta ^{})+2(\mu ^{\prime \prime }2\varphi ^{\prime \prime }),$$

$$4\varphi ^{\mathrm{\hspace{0.17em}2}}+[2(1+ϵ)\varphi ^{}2\varphi ^{}][2(1+ϵ)\varphi ^{}2(1+2ϵ)\varphi ^{}]+2[2(1+ϵ)\varphi ^{\prime \prime }2\varphi ^{\prime \prime }],$$

$$4\varphi ^{\mathrm{\hspace{0.17em}2}}+[2ϵ\varphi ^{}][2ϵ\varphi ^{}]+4\varphi ^{\prime \prime },$$

$$4\varphi ^{2}(1ϵ^2)+4ϵ\varphi ^{\prime \prime }$$

Left-hand side of Equation (8):

$$\mu ^{}(\mu ^{}4\varphi ^{}),$$

$$2(1+ϵ)\varphi ^{}[2(1+ϵ)\varphi ^{}4\varphi ^{}],$$

$$2(1+ϵ)\varphi ^{}[2\varphi ^{}+2ϵ\varphi ^{}],$$

$$4\varphi ^{2}(1ϵ^2).$$

The factor $`(1ϵ^2)`$ appears in the left-hand side of each equation multiplying $`\varphi ^{\mathrm{\hspace{0.17em}2}}`$. To satisfy Equation (4), one must have $`(1ϵ^2)=0`$, leading to the values $`ϵ=\pm 1`$ for the general relativity solutions. Equations (9) and (10) in the text are obtained after dividing the left-hand sides above by 8.

\[Note 3\] The only possible demand Fackerell can have would be, that the matter tensor should be $`\tau _\mu ^\nu =e^\eta \mathrm{diag}[(\sigma ,0,0,0)]`$ for Newtonian correspondence. Ah $`\mathrm{}`$, but that requires $`ϵ=0`$, which is precisely the Yılmaz Theory! Some general relativists refer to the matter tensor $`\tau _\mu ^\nu =e^\eta \mathrm{diag}[(1+ϵ,ϵ/2,ϵ/2,0)\sigma ]`$, with $`ϵ=\pm 1`$, shown here to be a solution of the general relativity field equations for slabs, as “peculiar matter.” However, the standard non-isotropic Schwarzschild metric solutions of general relativity have matter tensors of exactly this form!

\[Note 4\] To see that the Taub metric is devoid of physical content, substitute into the left-hand sides of the field equations (6) - (8) the relations which hold in the Taub metric: $`\mu =4\varphi `$, $`\eta =2\varphi `$, $`\varphi ^{\prime \prime }=4\varphi ^{\mathrm{\hspace{0.17em}2}}`$, to obtain:
Left-hand side of Equation (6):

$$\mu ^{}(3\mu ^{}2\eta ^{})+4\mu ^{\prime \prime },$$

$$4\varphi ^{}(12\varphi ^{}+4\varphi ^{})+16\varphi ^{\prime \prime },$$

$$64\varphi ^{2}64\varphi ^{2}0.$$

Left-hand side of Equation (7):

$$4\varphi ^{\mathrm{\hspace{0.17em}2}}+(\mu ^{}2\varphi ^{})(\mu ^{}\eta ^{})+2(\mu ^{\prime \prime }2\varphi ^{\prime \prime }),$$

$$4\varphi ^{\mathrm{\hspace{0.17em}2}}+(4\varphi ^{}2\varphi ^{})(4\varphi ^{}+2\varphi ^{})+2(4\varphi ^{\prime \prime }2\varphi ^{\prime \prime }),$$

$$4\varphi ^{\mathrm{\hspace{0.17em}2}}+12\varphi ^{2}+4\varphi ^{\prime \prime },$$

$$16\varphi ^{\mathrm{\hspace{0.17em}2}}16\varphi ^{\mathrm{\hspace{0.17em}2}}0.$$

Left-hand side of Equation (8):

$$\mu ^{}(\mu ^{}4\varphi ^{}),$$

$$4\varphi ^{}(4\varphi ^{}4\varphi ^{})0.$$

\[Note 5\] The issue of interactive N-body solutions is discussed quite generally by C.O. Alley, D. Leiter, Y. Mizobuchi and H. Yılmaz, “Energy Crisis in Astrophysics (Black Holes vs. N-Body Metrics),” available at the Los Alamos e-print archive, xxx.lanl.gov, as astro-ph/9906458 (28 June 1999). The absence of interaction between the two slabs is a specific example, admitting a simple analytic solution, of the situation in general relativity.