Phonons in Random Elastic Media and the Boson Peak
We show that the density of states of random wave equations, normalized by the square of the frequency, has a peak --- sometimes narrow and sometimes broad --- in the range of wave vectors between the disorder correlation length and the interatomic spacing. The results of this letter may be relevant for understanding vibrational spectra and light propagation in disordered solids.
peak10.tex
%% ****** Start of file template.aps ****** %
%%
%%
%%   This file is part of the APS files in the REVTeX 4 distribution.
%%   Version 4.0 of REVTeX, August 2001
%%
%%
%%   Copyright (c) 2001 The American Physical Society.
%%
%%   See the REVTeX 4 README file for restrictions and more information.
%%
%
% This is a template for producing manuscripts for use with REVTEX 4.0
% Copy this file to another name and then work on that file.
% That way, you always have this original template file to use.
%
% Group addresses by affiliation; use superscriptaddress for long
% author lists, or if there are many overlapping affiliations.
% For Phys. Rev. appearance, change preprint to twocolumn.
% Choose pra, prb, prc, prd, pre, prl, prstab, or rmp for journal
%  Add 'draft' option to mark overfull boxes with black boxes
%  Add 'showpacs' option to make PACS codes appear
%  Add 'showkeys' option to make keywords appear
%\documentclass[aps,prl,preprint,groupedaddress]{revtex4}
%\documentclass[aps,prl,preprint,superscriptaddress]{revtex4}
\documentclass[aps,prl,twocolumn,groupedaddress,showpacs]{revtex4}
\usepackage{epsfig}
\usepackage{graphicx}
\usepackage{amsfonts}
\usepackage{latexsym}
\usepackage{bm}
%\usepackage{pdfsync}


% You should use BibTeX and apsrev.bst for references
% Choosing a journal automatically selects the correct APS
% BibTeX style file (bst file), so only uncomment the line
% below if necessary.
%\bibliographystyle{apsrev}

\def\kx{\left( {p k_x \over 2} \right)}
\def\ky{\left( {p k_y \over 2} \right)}
\def\bk{{\bf k}}
\def\bx{{\bf x}}
\def\by{{\bf y}}
\newcommand{\boldgreek}[1]{\mbox{\boldmath$\mathbf{#1}$\unboldmath}}
\def\bpartial{\boldgreek{\partial}}
\def\8{\infty}
\def\oh{\frac{1}{2}}
\def\ot{\frac{1}{3}}
\def\oq{\frac{1}{4}}
\def\tt{\frac{2}{3}}
\def\ft{\frac{4}{3}}
\def\tq{\frac{3}{4}}
\def\d{\partial}
\def\i{\imath\,}
\def\ih{\frac{\imath}{2}\,}
\def\undertext#1{\vtop{\hbox{#1}\kern 1pt \hrule}}
\def\ra{\rightarrow}
\def\lfa{\leftarrow}
\def\Ra{\Rightarrow}
\def\lra{\longrightarrow}
\def\ler{\leftrightarrow}
\def\lrb#1{\left(#1\right)}
\def\O#1{O\left(#1\right)}
\def\VEV#1{\left\langle\,#1\,\right\rangle}
\def\tr{\hbox{tr}\,}
\def\trb#1{\tr\lrb{#1}}
\def\dd#1{\frac{d}{d#1}}
\def\dbyd#1#2{\frac{d#1}{d#2}}
\def\pp#1{\frac{\partial}{\partial#1}}
\def\pbyp#1#2{\frac{\partial#1}{\partial#2}}
\def\ff#1{\frac{\delta}{\delta#1}}
\def\fbyf#1#2{\frac{\delta#1}{\delta#2}}
\def\pd#1{\partial_{#1}}
\def\br{\\ \nonumber & &}
\def\brr{\right. \\ \nonumber & &\left.}
\def\inv#1{\frac{1}{#1}}
\def\be{\begin{equation}}
\def\ee{\end{equation}}
\def\bea{\begin{eqnarray} & &}
\def\eea{\end{eqnarray}}
\def\ct#1{\cite{#1}}
\def\rf#1{(\ref{#1})}
\def\EXP#1{\exp\left(#1\right)}
\def\TEXP#1{\hat{T}\exp\left(#1\right)}
\def\INT#1#2{\int_{#1}^{#2}}
\def\MAT{{\it Mathematica }}
\def\LHS{left-hand side }
\def\RHS{right-hand side }
\def\COM#1#2{\left\lbrack #1\,,\,#2\right\rbrack}
\def\AC#1#2{\left\lbrace #1\,,\,#2\right\rbrace}


\def\cH{{\cal H}}
\def\rf#1{(\ref{#1})}
\def\t{\tilde}
\def\cE{{\cal E}}

\def\rfs#1{Eq.~\rf{#1}}

\begin{document}

% Use the \preprint command to place your local institutional report
% number in the upper righthand corner of the title page in preprint mode.
% Multiple \preprint commands are allowed.
% Use the 'preprintnumbers' class option to override journal defaults
% to display numbers if necessary
%\preprint{}

%Title of paper
\title{Phonons in Random Elastic Media and the Boson Peak}

% repeat the \author .. \affiliation  etc. as needed
% \email, \thanks, \homepage, \altaffiliation all apply to the current
% author. Explanatory text should go in the []'s, actual e-mail
% address or url should go in the {}'s for \email and \homepage.
% Please use the appropriate macro foreach each type of information

% \affiliation command applies to all authors since the last
% \affiliation command. The \affiliation command should follow the
% other information
% \affiliation can be followed by \email, \homepage, \thanks as well.
\author{V. Gurarie}
%\email[]{Your e-mail address}
%\homepage[]{Your web page}
%\thanks{}
%\altaffiliation{}
\affiliation{Department of Physics, CB390, University of Colorado,
Boulder CO 80309, USA}
\author{A. Altland}
\affiliation{Institut f\"ur Theoretische Physik, Universit\"at zu
K\"oln, Z\"ulpicher Str 77, K\"oln, Germany}

%Collaboration name if desired (requires use of superscriptaddress
%option in \documentclass). \noaffiliation is required (may also be
%used with the \author command).
%\collaboration can be followed by \email, \homepage, \thanks as well.
%\collaboration{}
%\noaffiliation

\date{\today}

\begin{abstract}
We show that the density of states of random wave equations,
normalized by the square of the frequency,  has a peak ---
sometimes narrow and sometimes broad --- in the range of wave
vectors between the disorder correlation length and the
interatomic spacing. The results of this letter may be relevant
for understanding vibrational spectra and light propagation in
disordered solids.
\end{abstract}

% insert suggested PACS numbers in braces on next line
\pacs{63.50.+x}
% insert suggested keywords - APS authors don't need to do this
%\keywords{}

%\maketitle must follow title, authors, abstract, \pacs, and \keywords
\maketitle

% body of paper here - Use proper section commands
One of the intriguing features of random elastic media ---
observable in both Raman or neutron cross sections but also in
calorimetric measurements --- is an anomalous accumulation of
phonons at low frequencies \cite{Jackle1981}. This phenomenon
finds its most prominent manifestation in a peak in the quantity
$I(\omega)=\rho(\omega)/\omega^2$, where $\rho(\omega)$ is the
phonon density of states (DoS) and $\omega$ is the frequency. In
recent years, many competing theories as to the origin of this
so--called boson peak have been formulated (for a recent
reference, see Ref.~\cite{Parshin2003} and references therein). Of
these approaches a majority is based on model mechanisms specific
to the low temperature physics of amorphous solids. Yet the boson
peak is shown by both glassy and random crystalline systems alike,
an observation which has ignited the search for an explanation
which is not tied to the specifics of a glassy environment.
%\cite{Schirmacher2004}.

On a basic level, the acoustic excitations of both amorphous
materials and disordered crystals are described by random wave
equations. Existing analyses of such equations in the literature
indeed predicted a disorder generated excess DoS,
\cite{Schirmacher1998}. However, these structures where observed
at high frequencies (wavelengths of the order of the interatomic
spacing), while the boson peak is a low energy phenomenon. In this
letter we argue that the DoS of elastic vibrations in disordered
media is enhanced by disorder at low frequencies, with $I(\omega)$
exhibiting either a peak or a broad maximum (see Fig.~\ref{graph1}
for a representative picture of the DoS). Which of the two will be
observed crucially depends on the modelling of the randomness, i.
e. the \textit{type} of disorder at work. Since the profile of the
randomness of `real' systems is generally unknown, we are not in a
position to judge whether the present mechanism alone may account
for the spectral peaks observed in all experiments. We believe,
however, that it operates under quite general conditions.

%% \begin{figure}[htbp]
%% \centerline{\epsfxsize=3.5in  \epsfbox{dos.eps} }
%% \caption{\label{graph1} $I(\omega)$ for the random elastic media
%% in cases where $\alpha_2>0$ or $\alpha_2<0$.}
%% \end{figure}

\begin{figure}[htbp]
\centerline{\resizebox{3.5in}{!}{\includegraphics{dos}} }
\caption{\label{graph1} $I(\omega)$ for the random elastic media
in cases where $\alpha_2>0$ or $\alpha_2<0$.}
\end{figure}

Let us begin by discussing some generic large scale structures of
the DoS of acoustic excitations. In a clean system, the linear
relation $\omega = \bar c k$ between frequency and magnitude $k$
of the wave vector gives rise to a DoS
\begin{equation}
\label{eq:naiveDoS} \rho(\omega) = \int {d^d k \over (2 \pi)^d
}~\delta(\omega-\bar c k) = {A_d \over (2 \pi \bar c)^d }~
\omega^{d-1}.
\end{equation}
Here, $A_d$ is the area of the $d$--dimensional unit sphere and
$\bar c$ the speed of sound (for our present analysis, the
potential existence of several acoustic phonon branches with
different sound velocities will not be of importance). Tied to the
linearizability of the dispersion relation, the relation $\rho
\sim \omega^{d-1}$ applies only to  frequencies $\omega\ll
\omega_D$ much smaller than the Debye frequency $\omega_D\sim
{\bar c\over a}$ ($a$: lattice spacing.) On dimensional grounds,
the generalization of (\ref{eq:naiveDoS}) to higher frequencies
must be of the form
\begin{equation}
\label{eq:latt} \rho(\omega) \sim \omega^{d-1} \tilde
f\left(\omega/
    \omega_D\right),
\end{equation}
where $\tilde f(u)$ is some function with $\tilde f(u\ll 1)\simeq
1$. It is clear  that $\t f(u)$ falls off to zero at $u \gg 1$.

In contrast, in disordered materials we find strong deviations
from the above scaling form already for frequencies $\omega\ll
\omega_D$. Indeed, our analysis of the random wave equation below
will lead us to
\begin{equation}
  \label{eq:DoSscalling}
  \rho(\omega)= {A_d \over (2 \pi \bar c)^d }~\omega^{d-1}
  f\left({\omega l\over \bar c}\right),\;\;
  \omega \ll \omega_D,
\end{equation}
where the effective correlation length $l\gg a$ of the disorder is
assumed to be much larger than the lattice spacing, and $\bar c$ is
the typical speed of sound in random media to be defined more
precisely below. The scaling function $f$ is given by
\begin{eqnarray}
\label{eq:fexp}
f(u) &=& 1 + \alpha_1 u^2 + \dots, \ u \ll 1, \cr
f(u) &=& \beta \left( 1 + {\alpha_2 \over u^2} + \dots \right),
\;1\ll  u \ll {l \over a},
\end{eqnarray}
where the form of the coefficients $\alpha_{1,2}$, $\beta$ depends
on the space dimensionality and on a few basic characteristics of
the randomness. In the most interesting case, $d=3$, we find that
$\beta>1$ and $\alpha_1>0$. The sign of $\alpha_2$ lacks
universality. $\alpha_2$ is negative if disorder in materials is
mostly in the random elastic constants, and it is positive if
disorder is mostly due to fluctuating mass density.

%Using a discrete language --- wherein one should
%think of a system of `random masses coupled by random spring
%constants' ---, $\alpha_2>0$ ($\alpha_2<0$) if fluctuations in the
%mass (the spring constants) are the dominant source of randomness.
For $\alpha_2>0$, $I(\omega)$ grows for $\omega \ll \bar c/l$ and
falls off at $\omega \gg \bar c/l$ (while still $\omega \ll
\omega_D$). This implies the existence a peak
at $\omega \sim \bar c/l$. In contrast, for $\alpha_2<0$,
$I(\omega)$ increases until $\omega$ becomes much larger than
$\bar c/l$.  Combined with the expected drop--off of $I(\omega)$
at $\omega \gtrsim \omega_D$, this produces a broad maximum for
$I(\omega)$ between $\bar c/l$ and $\omega_D$. Fig.~\ref{graph1}
shows a caricature of the two scenarios. We finally
note that if the disorder is relatively weak, ${\rm rms\,}c<\bar
c$, (so that the fluctuations of local speed of sound $c$ are
smaller than the typical sound velocity) $\alpha_1\sim
\alpha_2\sim {\,\rm var}\, c/\bar c^2$; for more generally
applicable expressions, see below.



%Referring for a discussion of other space dimensionalities to [],
%we here merely note that for $d=1$ and certain types of the
%disorder distribution the DoS can be computed by exact
%methods~\cite{Gurarie2002a}. The results thus obtained conform
%with  the profile of our scaling function above.


%% Here at $d=3$, $\beta>1$ and $\alpha_1>0$, $\alpha_2 >0$.  When $d
%% \not = 3$, the behavior of $f(u)$ is quite different, however.  At
%% $d=1$, $\beta<1$, $\alpha_1<0$ and $\alpha_2>0$, which is consistent
%% with the density of states being a monotonously decreasing function.
%% This can be confirmed by the exact density of states which can be
%% found at $d=1$ for some disorder distributions \cite{Gurarie2002a}. At
%% $d=2$, $\beta=1$, $\alpha_2=0$ and there is evidence, which we did not
%% explore in this letter, that the expansion \rfs{eq:fexp} acquires
%% logarithmic corrections. If $d>3$, then $\beta>1$ and $\alpha_1>0$.
%% However, $\alpha_2=0$ at $d=4$ and it becomes negative at $d>4$.

We next turn to the derivation of these results. Consider the
wave equation
\begin{equation}
  \label{eq:random_wave}
  \left(\Delta + {\omega^2 \over c^2({\bf x})}\right) \psi({\bf x})=0,
\end{equation}
where the random velocity field $c({\bf x})$ fluctuates on spatial
scales $\sim l$.  Interpreting the variable $m({\bf x})\equiv
1/c^2({\bf x})$ as the mass density of the random medium, we refer to
\rfs{eq:random_wave} as a `random mass' wave equation.

In the regime of small frequencies, $\omega \ll \bar c/l$, we are
facing a situation where the correlation length $l$ is much smaller
than the `typical' wave length $k^{-1} \equiv \bar c/\omega$ at which
the wave function $\psi({\bf x})$ fluctuates.  The phonon field
effectively averages over many fluctuation intervals of the disorder
and, to a first approximation, \rfs{eq:random_wave} may be replaced
for its average over random $c({\bf x})$: $\left(\Delta + {\omega^2
\over \bar c^2}\right) \psi({\bf x})=0$, where the typical speed of
sound $\bar c$ is defined as
\begin{equation}
\label{eq:typc} \bar c \equiv \langle c({\bf
x})^{-2}\rangle^{-1/2}
\end{equation}
and the angular brackets denote averaging over random $c({\bf x})$.
The very low frequency asymptotics is then given by \rfs{eq:naiveDoS}
with $\bar c$ taken from \rfs{eq:typc}.  (This result, as well as the
high frequency asymptotics \rfs{eq:beta} discussed below, are
originally due to Chalker \cite{ChalkerPrivate}.)


 To compute corrections to the low frequency asymptotics of the
DoS we rewrite Eq.~(\ref{eq:random_wave}) as
\begin{equation}
  \label{eq:rw_pert}
    \left(- \Delta - {\omega^2\over \bar c^2} h(\bx)\right) \psi({\bf x})= {\omega^2
    \over \bar c^2} \psi(\bx),
\end{equation}
where the function $ h({\bf x})\equiv \left[{\bar c/ c({\bf
x})}\right]^2-1 $ describes the randomness. This representation
suggests to interpret $\omega^2/\bar c^2$ as the eigenvalue of the
operator $\Delta$ weakly perturbed by $\omega^2 h(\bx)/\bar c^2$.
%(Two remarks: (i) To lowest order perturbation theory, the factor
%$\omega^2$ multiplying $h$ can be replaced by the unperturbed
%eigenvalue. (ii) In a sense to be clarified below, the `weakness'
%of the perturbation follows from the fact that the random variable
%$h$ has zero average.)
The unperturbed problem is trivially
diagonalized by a set of plane waves $ \psi_\bk(\bx) \equiv
{1\over L^{d/2}} \exp(i\bk \cdot \bx), $ with eigenvalues
$\omega^2/\bar c^2 = k^2$ ($L$ is the system size.) We next apply
standard perturbation theory to compute the eigenvalue shift
caused by the presence of the perturbation $\omega^2 h/\bar c^2$.
The vanishing of $\VEV{h}$ implies that, on average, there are no
first order corrections. To second order in $h$ we find that the
average eigenvalue is given by %%
 % $$
 % {\omega^2\over \bar c^2} = k^2 + {\omega^4\over \bar c^4 L^{2d}}
 % \sum_{\bk'} {g(\bk-\bk')\over
 % k^2-k'^2},
 % $$
 %%
$
  {\omega^2\over \bar c^2} = k^2 + k^4 \int {d^d
k' \over (2 \pi)^d} {g({\bf k'}) \over k^2-k'^2}$,
where $g({\bf k})$ is the Fourier transform of the disorder
correlation function of $g({\bf x}-{\bf
y})\equiv \VEV{h({\bf x}) h({\bf y})}$. Irrespective of the
distribution of the disorder, this function (a) drops off rapidly for
$k>l^{-1}$ and (b) approaches a constant value $\sim  l^d {\rm
\,var}(h)$ for $k\ll l^{-1}$. This implies that
for $d> 2$, the integral above is dominated by momenta $k'\sim l^{-1}$.
At the same time, the reference momentum $k\sim \omega/\bar c\ll
l^{-1}$ is small. We thus neglect the $k$--dependence of the
integrand and arrive at
\begin{equation}
  \label{eq:eval_lowo}
  {\omega^2\over \bar c^2} = k^2 -  k^4 \int {d^d
k' \over (2 \pi)^d} {g({\bf k'}) \over k'^2}.
\end{equation}
%At this point a remark is in order:
Physically, the correction to the
zeroth order eigenvalue is due to virtual scattering events wherein
states with low--lying momentum $k\sim \omega/\bar c$ scatter off
rapid fluctuations of $h$ into high--lying states $k'\sim l^{-1}$.
Yet for small frequencies, the large phase volume $\sim l^{-d}$
available to these scattering processes cannot outweigh the overall
multiplicative factor $\omega^4$.  This mechanism pervades to higher
orders in the eigenvalue expansion and justifies the perturbative
approach.  In particular, Eq.  (\ref{eq:eval_lowo}) indeed describes
the dominant correction to the low frequency dispersion relation.

We
next substitute Eq.  (\ref{eq:eval_lowo}) into $\rho(\omega) = {1\over
L^d} \sum_\bk \delta(\omega(k)-\omega)$ to obtain an expansion of the
DoS as in Eqs.  (\ref{eq:DoSscalling},\ref{eq:fexp}).  Specifically,
the coefficient
\begin{equation} \label{eq:alpha1} \alpha_1 ={d+2
\over 2l^2} \int {d^d k \over (2 \pi)^d} {g({\bf k}) \over k^2}\sim
\mathrm{var}(h)>0,
\end{equation}
where the proportionality to $\mathrm{var}(h)$ follows from the
properties of the correlation function $g$ discussed above.
Interestingly, positivity of $\alpha_1$ in \rfs{eq:alpha1} is a
direct consequence of the fact that the second order perturbation
theory always lowers the ground state energy.

For small velocity fluctuations, ${\rm \,var}(c) \ll \bar c^2$, we
obtain the estimate $\alpha_1\sim \mathrm{var}(h)\sim{\rm
var}(c)/\bar c^2$ quoted above. We finally note that these results
hold only for $d>2$. For $d \le 2$, the integral in
\rfs{eq:alpha1} is infrared divergent and the approximation scheme
employed here breaks down.

% We would also like to note that there exists an alternative way of
% deriving \rfs{eq:alpha1} which allows for a systematic expansion
% of the DoS in powers of $\omega$. Introduce a Green's function
% $$
% G = {1 \over \bar c^2} \left[ -\Delta^2 - {\omega^2 \over \bar
% c^2}\left(1+h({\bf x}) \right)- i \epsilon \right]^{-1}
% \left(1+h({\bf x}) \right).
% $$
% The DoS is related to it in a trivial fashion, $\rho(\omega) =(2
% \omega /\pi)~{\rm Im}\VEV{G(x,x)}$. Expanding $G$ in powers of $h$
% and averaging over disorder term by term, it is possible, in
% principle, to calculate the low frequency expansion of the density
% of states up to arbitrary order. In the lowest order this method,
% of course, reproduces \rfs{eq:alpha1}.


We next turn to the discussion of the high frequency case, $\omega
\gg \bar c/l$. In this regime, the velocity field varies very
little over length scales comparable to the typical wavelength.
This implies that, locally, the solutions of
(\ref{eq:random_wave}) behave like plane waves with the local
dispersion relation $k({\bf x})=\omega/c({\bf x})$ and the local
density of states  given by
$$
\rho(\omega,{\bf x}) = {A_d \over (2 \pi)^d c^d({\bf x})}~
\omega^{d-1}.
$$
Averaging this result over the random velocity field we obtain the
global density of states
$$
\rho(\omega) = {A_d \over (2 \pi)^d} \VEV{1 \over c^d({\bf
x})}~\omega^{d-1}, \qquad\omega \gg {\bar c \over l}.
$$
We thus find that the coefficient $\beta$, introduced in
\rfs{eq:fexp}, is given by
\begin{equation}
\label{eq:beta} \beta= {\VEV{c^{-d} \left( {\bf x} \right)} \over
\VEV{c^{-2}({\bf x})}^{d \over 2}}.
\end{equation}
As a consequence of the convexity of the power function, $\beta>1$
for $d>2$.


 To obtain corrections to the high frequency
 asymptotics~(\ref{eq:beta}), we
 need
to compute distortions in the spectral density caused by shallow
(compared to the wave length) variations  of the velocity field.
It turns out that this task is most efficiently tackled by
analyzing the Green's function
\begin{equation}
\label{eq:greens} G=-\left( {1\over \omega^2} \nabla^2 + m(\bx) +
i \epsilon \right)^{-1},
\end{equation}
where $m\equiv c^{-2}$. From (\ref{eq:greens}) the
average DoS is obtained as
\begin{equation}
\label{eq:dosgr} \rho(\omega) =  {2  \over \pi \omega} ~{\rm Im}
\left\langle m({\bf x})~G(\bx,\bx)\right\rangle,
\end{equation}
What makes the operator (\ref{eq:greens}) a good  starting point
for our analysis is its structural similarity to a Schr\"odinger
operator with `Planck's constant' $\hbar\sim \omega^{-1}$. This
analogy will enable us to apply semiclassical approximation
schemes familiar from quantum mechanics.  We begin by applying the
Wigner transform
$$
G(\bx,\by) = {\omega^d \over (2 \pi)^d} \int d^d k~G(\bx,\bk)~e^{i
\omega
  \bk\cdot \left(\bx-\by\right)},
$$
whereupon the `Schr\"odinger equation' assumes the form
\begin{equation}
\label{eq:Schro} \left[\left({\bf k}-i\omega^{-1}\partial_{\bf
x}\right)^2 -m({\bf x})-i\epsilon\right]G(\bx,\bk)=1.
\end{equation}
To make use of the smallness, $\sim \omega^{-1}$, of the
derivative operators, we expand $G(\bx,\bk)$ in powers of
$\omega^{-1}$,
$$
G({\bf x},{\bf k}) = G^{(0)}({\bf x},{\bf k})+ \frac{1}{\omega}
G^{(1)}({\bf x},{\bf k}) + \frac{1}{\omega^2} G^{(2)}({\bf x},{\bf
k})+ \dots,
$$
substitute the expansion into \rfs{eq:Schro}, and find $G^{(0)}$,
$G^{(1)}$, $G^{(2)}$, $\dots$, recursively. By symmetry,
$G^{(1)}$ vanishes upon configurational averaging so that the
dominant correction to the DoS is provided by $G^{(2)}$.
Substituting this coefficient into \rfs{eq:dosgr}, and comparing
with Eq.~(\ref{eq:fexp}) we obtain the result
$$
\alpha_2 ={\left(d-2\right)^2 \left(4-d \right) \over 24} { l^2
\over \bar c^2 \VEV{c^{-d}}}  \VEV{{ \nabla c \cdot \nabla c \over
c^d}}.
$$
Clearly, $\alpha_2>0$ if $d=3$. If the fluctuations of $c$ are weak,
$\alpha_2 \sim {\rm var}\,c/\bar
c^2\sim \alpha_1$.

We thus derived \rfs{eq:fexp}. The positivity of the two expansion
coefficients $\alpha_{1,2}$ implies that the scaling function
$f(u)$ in \rfs{eq:DoSscalling} grows at $u \ll 1$ and falls off at
$u \gg 1$. This implies that the functional profile of the DoS
contains a peak somewhere at $u =\mathcal{O}(1)$.

To conclude our analysis of the prototypical wave equation
(\ref{eq:random_wave}), let us briefly discuss the case of space
dimensions $d\not=3$. For $d=1$, the high frequency coefficient
$\alpha_2>0$, yet $\beta<1$. Due to the instability in the the low
frequency expansion, the methods applied in this letter do not
allow to calculate $\alpha_1$ at $d \le 2$. However, for certain
distributions of the disorder transfer matrix methods may be
applied to obtain an exact solution~\cite{Gurarie2002a}. These
calculations show that $f(u)$ is a globally decreasing function,
and $\alpha_1<0$. For $d=2$, the two dominant high--frequency
coefficients are structureless, $\beta=1$ and $\alpha_2=0$, while
the low frequency expansion in \rfs{eq:fexp} will now involve
terms proportional to $\log(u)$. From these results we cannot
decide whether $\rho(\omega)/\omega$ has a maximum or a minimum
between its low and high frequency asymptotics. Finally, for
$d>3$, $\alpha_1>0$ and $\beta>1$, however $\alpha_2\le 0$, which
implies that $f(u)$ is a monotonously increasing function.
Summarizing, we see that  the normalized DoS of
\rfs{eq:random_wave}, $\rho(\omega)/\omega^{d-1}$, exhibits a peak
at wavelengths of the order of disorder correlation length only
for $d=3$.


Before proceeding, let us briefly compare our so far results to
earlier work.  In most numerical simulations of the problem (cf.
e.g., Ref.~\cite{Schirmacher1998} and references therein), the
disorder is chosen to be uncorrelated, that is, its correlation
length is of the order of the lattice spacing.  For such type of
disorder, $u\simeq 1$ translates to $\omega\sim \omega_D$ deep in
the bulk of the spectrum.  At these frequencies it is hard to tell
whether deviations from the low frequency asymptotics
$\rho(\omega)\sim \omega^2$ are caused by lattice effects (cf.
Eq.~(\ref{eq:latt})) or by disorder (Eq.~(\ref{eq:DoSscalling})).
Qualitatively, however, the numerical data is in agreement with
the results of our present analysis.  Turning to earlier
analytical work, we notice that most approaches to random elastic
problems rely on the self consistent Born approximation (SCBA).
(See Ref.~\cite{EfetovBook} for a general review of the methods
involved and Refs.~\cite{John1982,Schirmacher2002}  for the
applications of these methods to random wave equations.)  However,
for a number of reasons the results obtained by this method lack
quantitative reliability when used to calculate the DoS: (i)
Within the standard SCBA, $1/c^2$ is modelled as a Gaussian
distributed variable. This implies the existence of rare domains
where $1/c^2$ is negative. Using the techniques introduced in
Ref.~\cite{Gurarie2003b} it is possible to show that these
domains, no matter how small, lead to unstable wave modes and to
unphysical sharply growing contributions to the density of states
at higher frequencies.
%\footnote{The problem of unstable
%domains can be circumvented by treating $1/c$ rather than $1/c^2$ as
%Gaussian distributed.  This technique proves useful in analyzing
%\rfs{eq:bond1}, discussed below.}
(ii) The SCBA neglects certain contributions (`crossed diagrams')
to the perturbative expansion of the Green function
(\ref{eq:greens}).  At high frequencies, this approximation
becomes invalid. The value of the threshold frequency beyond which
the SCBA breaks down depends on the type of disorder under
consideration. At any rate, however, it is parametrically smaller
than the band width. At the same time, (iii), previous studies
focused on the case of short range correlated disorder. In this
case, deviations of the scaling law $\rho\sim \omega^{d-1}$ are
expected at frequencies $\omega \sim \omega_D$, i.e. in a regime
where the SCBA breaks down {\it and} lattice effects intervene. (The alternative techniques in
this letter avoid all of these problems.)



The wave equation \rfs{eq:random_wave} applies to the specific
case, where only the mass density of the elastic medium
fluctuates. The continuum description of a more general
environment wherein the elastic constants also fluctuate reads as
\begin{equation}
\label{eq:bond1} \left( \nabla \mu({\bf x}) \nabla + \omega^2
m({\bf x}) \right) \psi({\bf x}) =0.
\end{equation}
Here both $\mu({\bf x})$ and $m({\bf x})$ are random positive
quantities.

It turns out that the direct perturbative expansion applied above
to \rfs{eq:random_wave} cannot be used to determine the
coefficient $\alpha_1$ of the problem \rfs{eq:bond1}. As an
alternative, we apply a generalized variant of the
self--consistent Born approximation, wherein
$\mu(\bx)=\mu_0+\sigma^2(\bx)$, and $\sigma(\bx)$ is a Gaussian
distributed variable with zero mean. (In this way, positivity of
the elastic constant is ensured.) The actual implementation of the
SCBA for this type of disorder turns out to be rather involved and
its details will be published elsewhere. Suffice to say that at
$d=3$, $\alpha_1$ is still positive. $\alpha_2$, on the other
hand, no longer has definite sign. Applying the high frequency
expansion outlined above for the random mass density case we find
\begin{eqnarray}
\label{eq:bond} \alpha_2 ={l^2 \over \bar c^2 \VEV{m^{3 \over 2}
\mu^{-{3 \over 2}}}} \left[{1 \over 96} \VEV{ {\mu^2 \left( \nabla
m \right)^2
\over  m^{3 \over 2} \mu^{5 \over 2}}} + \right. \\
\nonumber \left. {5 \over 48} \VEV{m \mu \nabla m \nabla \mu \over
m^{3 \over 2} \mu^{5 \over 2}} -{23 \over 96}\VEV {m^2
\left(\nabla \mu \right)^2 \over m^{3 \over 2} \mu^{5 \over
2}}\right].
\end{eqnarray}
If $\alpha_2$ is negative (which happens, for example, in the
limiting case of non--random $m({\bf x})$) $f(u)$ no longer has a
maximum at $u \simeq 1$. Instead, $I(\omega)$ is an increasing
function of $\omega$ with a broad maximum reached at frequencies
$\omega > \bar c/l$, before dropping off at frequencies higher
than $\omega_D$. Concluding, we find that, depending on the
profile of the disorder,
 the envelope function $I(\omega)$ of the random mass/elastic constants problem
 either contains a low--frequency peak, or a broad high frequency
 maximum.


So far, we have considered the case of scalar phonons.
%$\psi({\bf
%  x})$.
In a realistic environment, however, $\psi\to u_i$ will be a
  $d$--component {\it vector}. The most general random vector problem
  would be governed by a formidable rank four random elastic modulus
  tensor. Assuming, however, a medium consisting of a
  random accumulation of `micro--crystallites' each of which
  possessing intrinsic rotational invariance --- this assumption may
  well be violated, especially in glassy environments, although it seems
  to work in polycrystalline materials --- the
  effective wave equation reduces to
\begin{equation}
\label{eq:elastic} \d_i \left[\lambda({\bf x}) \d_j u_j\right] +
\d_j \left[ \mu({\bf x}) \left( \d_i u_j + \d_j u_i \right)\right]
+ \omega^2 m({\bf x}) u_i =0,
\end{equation}
where $m(\bx)$ is a random density of the medium and $\lambda(\bx)$,
$\mu(\bx)$ are random Lam\'e coefficients.
%While the general analysis
%of Eq.(\ref{eq:elastic}) is still very complicated,
%It is now instructive to consider a few limiting cases.

In the limit where only $\lambda(\bx)$ is random, \rfs{eq:elastic}
can be mapped into \rfs{eq:random_wave} by substitution of
$\psi=\d_i u_i (\lambda + 2 \mu)$. In this case, both $\alpha_1$
and $\alpha_2$ are positive.
%We note that only
%longitudinal phonons interact with disorder, while the transverse
%ones propagate freely. \
This result carries over to the case where, in addition to random
$\lambda(\bx)$, also the density $m(\bx)$ is
made random; For these two types of disorder
$I(\omega)$ has a peak at $\omega \sim \bar c/l$.
If, however, the shear modulus $\mu$ is also random, \rfs{eq:elastic}
becomes more similar to \rfs{eq:bond1}. Applying the SCBA, it is still
possible to show that $\alpha_1>0$ (for $d>2$). As with our previous
discussion of (\ref{eq:bond1}), however, $\alpha_2$ is no longer of
definite sign. Specifically, $I(\omega)$ will contain a broad maximum
($\alpha_2<0$) if
 only the shear modulus $\mu(\bx)$ is random.


 To conclude, we have shown that, depending on the type of disorder,
 the normalized DoS $I(\omega)$ of random wave equations  may either contain a peak at phonon
 wavelengths of the order of the disorder correlation length, or a
 broad maximum at wavelengths below the correlation length. We believe
 that this work is not only relevant for the interpretation of data on
 vibrational modes in random media but also to the analysis
 of other types of waves, such as electromagnetic waves in random
 environments.

 The authors are grateful to J. T. Chalker for many valuable
 discussions. VG is also grateful to A. V. Andreev for useful
 comments. Work supported by SFB/TR 12 of the Deutsche Forschungsgemeinschaft.

%The high frequency asymptotics can be found if we note that in
%this regime $g(p)$ can effectively be replaced by $M (2 \pi)^3
%\delta(p)$, where $M$ sets the characteristic scale of the
%correlation function $g(p)$. In other words, all the length scales
%become much smaller than $l$, the disorder correlation length.
%Then $Q_{1,2} = M \lambda_{1,2}$ and we find
%\begin{equation}
%\lambda_1 = - {p^2 \over 2 M} - {i \over 2} \sqrt{  4 {p^2 \over
%M} \omega^2 - {p^4 \over M^2 }}.
%\end{equation}
%Therefore, the density of states is given by
%\begin{equation}
%\rho(\omega) = {1 \over \pi \omega^3} \int_0^{2 \omega \sqrt{M}}
%{p^2 dp \over 2 \pi^2} \sqrt{  4 {p^2 \over M} \omega^2 - {p^4
%\over M^2 }}
%\end{equation}







% Recently low frequency bosonic excitations in random media were
%discussed in a very general framework in
%Ref.~\cite{Gurarie2002,Gurarie2003a}. The discussion in these
%publications showed that the behavior of those excitations depend
%crucially on whether or not they are Goldstone modes. In turn, the
%Goldstone modes can be further subdivided into the phonon-like and
%magnon-like modes, depending on the structure of their
%Hamiltonians.

%The phonon-like modes in random elastic media were first
%thoroughly discussed two decades ago in Ref.~\cite{John1982}. One
%of the basic results of that publication was that in an elastic
%media with short range correlated random elastic parameters the
%low frequency density of states goes as
%\begin{equation}
%\label{eq:DoS} \rho(\omega) \propto \omega^{d-1}, \, \omega \ll
%\omega_{D},
%\end{equation}
%where $d$ is the dimensionality of the media, and $\omega_D$ is
%the Debye frequency. This result can be anticipated without any
%calculations. The long wavelength nature of low frequency phonons
%(and indeed of any Goldstone modes) tend to average local disorder
%over distances of the order of wavelength and therefore the
%fluctuations in the elastic properties of the media can be
%completely neglected. On the other hand, it is well known that the
%density of states in any clean elastic media follows \rfs{eq:DoS}.













% References should be done using the \cite, \ref, and \label commands
%\section{}
% Put \label in argument of \section for cross-referencing
%\section{\label{}}
%\subsection{}
%\subsubsection{}

% If in two-column mode, this environment will change to single-column
% format so that long equations can be displayed. Use
% sparingly.
%\begin{widetext}
% put long equation here
%\end{widetext}

% figures should be put into the text as floats.
% Use the graphics or graphicx packages (distributed with LaTeX2e)
% and the \includegraphics macro defined in those packages.
% See the LaTeX Graphics Companion by Michel Goosens, Sebastian Rahtz,
% and Frank Mittelbach for instance.
%
% Here is an example of the general form of a figure:
% Fill in the caption in the braces of the \caption{} command. Put the label
% that you will use with \ref{} command in the braces of the \label{} command.
% Use the figure* environment if the figure should span across the
% entire page. There is no need to do explicit centering.

% \begin{figure}
% \includegraphics{}%
% \caption{\label{}}
% \end{figure}

% Surround figure environment with turnpage environment for landscape
% figure
% \begin{turnpage}
% \begin{figure}
% \includegraphics{}%
% \caption{\label{}}
% \end{figure}
% \end{turnpage}

% tables should appear as floats within the text
%
% Here is an example of the general form of a table:
% Fill in the caption in the braces of the \caption{} command. Put the label
% that you will use with \ref{} command in the braces of the \label{} command.
% Insert the column specifiers (l, r, c, d, etc.) in the empty braces of the
% \begin{tabular}{} command.
% The ruledtabular enviroment adds doubled rules to table and sets a
% reasonable default table settings.
% Use the table* environment to get a full-width table in two-column
% Add \usepackage{longtable} and the longtable (or longtable*}
% environment for nicely formatted long tables. Or use the the [H]
% placement option to break a long table (with less control than
% in longtable).
% \begin{table}%[H] add [H] placement to break table across pages
% \caption{\label{}}
% \begin{ruledtabular}
% \begin{tabular}{}
% Lines of table here ending with \\
% \end{tabular}
% \end{ruledtabular}
% \end{table}

% Surround table environment with turnpage environment for landscape
% table
% \begin{turnpage}
% \begin{table}
% \caption{\label{}}
% \begin{ruledtabular}
% \begin{tabular}{}
% \end{tabular}
% \end{ruledtabular}
% \end{table}
% \end{turnpage}

% Specify following sections are appendices. Use \appendix* if there
% only one appendix.
%\appendix
%\section{}

% If you have acknowledgments, this puts in the proper section head.
%\begin{acknowledgments}
% put your acknowledgments here.
%\end{acknowledgments}

% Create the reference section using BibTeX:
\bibliography{ref}

\end{document}
%
% ****** End of file template.aps ******

peak10.bbl
\begin{thebibliography}{9}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname bibnamefont\endcsname\relax
  \def\bibnamefont#1{#1}\fi
\expandafter\ifx\csname bibfnamefont\endcsname\relax
  \def\bibfnamefont#1{#1}\fi
\expandafter\ifx\csname citenamefont\endcsname\relax
  \def\citenamefont#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\providecommand{\bibinfo}[2]{#2}
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem[{\citenamefont{J\"ackle}(1981)}]{Jackle1981}
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{J\"ackle}}, in
  \emph{\bibinfo{booktitle}{Amorphous Solids}}, edited by
  \bibinfo{editor}{\bibfnamefont{W.~A.} \bibnamefont{Phillips}}
  (\bibinfo{publisher}{Springer-Verlag}, \bibinfo{address}{Heidelberg},
  \bibinfo{year}{1981}), chap.~\bibinfo{chapter}{8}, p. \bibinfo{pages}{135}.

\bibitem[{\citenamefont{{Gurevich {\it et al.}}}(2003)}]{Parshin2003}
\bibinfo{author}{\bibfnamefont{V.~L.} \bibnamefont{{Gurevich {\it et al.}}}},
  \bibinfo{journal}{Phys. Rev. B} \textbf{\bibinfo{volume}{67}},
  \bibinfo{pages}{094203} (\bibinfo{year}{2003}).

\bibitem[{\citenamefont{{Schirmacher {\it et al.}}}(1998)}]{Schirmacher1998}
\bibinfo{author}{\bibfnamefont{W.}~\bibnamefont{{Schirmacher {\it et al.}}}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{81}},
  \bibinfo{pages}{136} (\bibinfo{year}{1998}).

\bibitem[{\citenamefont{Chalker}()}]{ChalkerPrivate}
\bibinfo{author}{\bibfnamefont{J.~T.} \bibnamefont{Chalker}}, \eprint{private
  communication}.

\bibitem[{\citenamefont{Gurarie}(2002)}]{Gurarie2002a}
\bibinfo{author}{\bibfnamefont{V.}~\bibnamefont{Gurarie}},
  \bibinfo{journal}{unpublished}  (\bibinfo{year}{2002}).

\bibitem[{\citenamefont{Efetov}(1997)}]{EfetovBook}
\bibinfo{author}{\bibfnamefont{K.~B.} \bibnamefont{Efetov}},
  \emph{\bibinfo{title}{Supersymmetry in Disorder and Chaos}}
  (\bibinfo{publisher}{Cambridge University Press},
  \bibinfo{address}{Cambridge, UK}, \bibinfo{year}{1997}).

\bibitem[{\citenamefont{{John {\it et al.}}}(1982)}]{John1982}
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{{John {\it et al.}}}},
  \bibinfo{journal}{Phys. Rev. B} \textbf{\bibinfo{volume}{27}},
  \bibinfo{pages}{5592} (\bibinfo{year}{1982}).

\bibitem[{\citenamefont{{Schirmacher {\it et al.}}}(2002)}]{Schirmacher2002}
\bibinfo{author}{\bibfnamefont{W.}~\bibnamefont{{Schirmacher {\it et al.}}}},
  \bibinfo{journal}{Phys. Stat. Sol. (b)} \textbf{\bibinfo{volume}{230}},
  \bibinfo{pages}{31} (\bibinfo{year}{2002}).

\bibitem[{\citenamefont{Gurarie and Altland}(2004)}]{Gurarie2003b}
\bibinfo{author}{\bibfnamefont{V.}~\bibnamefont{Gurarie}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Altland}},
  \bibinfo{journal}{J Phys A: Math. Gen.} \textbf{\bibinfo{volume}{37}},
  \bibinfo{pages}{9357} (\bibinfo{year}{2004}).

\end{thebibliography}
