# A single defect approximation for localized states on random lattices

## Abstract

Geometrical disorder is present in many physical situations giving rise to eigenvalue problems. The simplest case of diffusion on a random lattice with fluctuating site connectivities is studied analytically and by exact numerical diagonalizations. Localization of eigenmodes is shown to be induced by geometrical defects, that is sites with abnormally low or large connectivities. We expose a “single defect approximation” (SDA) scheme founded on this mechanism that provides an accurate quantitative description of both extended and localized regions of the spectrum. We then present a systematic diagrammatic expansion allowing to use SDA for finite-dimensional problems, e.g. to determine the localized harmonic modes of amorphous media.

Since Anderson’s fundamental work , physical systems in presence of disorder are well known to exhibit localization effects . While most attention has been paid so far to Hamiltonians with random potentials (e.g. stemming from impurities), there are situations in which disorder also originates from geometry.

Among these are of particular interest the harmonic vibrations of amorphous materials as liquids, colloids, glasses, … around random particle configurations. Recent experiments on sound propagation in granular media  have stressed the possible presence of localization effects, highly correlated with the microscopic structure of the sample. The existing theoretical framework for calculating the density of harmonic modes in amorphous systems was developed in liquid theory. In this context, microscopic configurations are not frozen but instantaneous normal modes (INM) give access to short time dynamics . Wu and Loring  and Wan and Stratt  have calculated good estimates of the density of INM for Lennard-Jones liquids, averaged over instantaneous particle configurations. However, localization-delocalization properties of the eigenvectors have not been considered.

Diffusion on random lattices is another problem where geometrical randomness plays a crucial role . Long time dynamics is deeply related to the small eigenvalues of the Laplacian on the lattice and therefore to its spectral dimension. Campbell suggested that diffusion on a random lattice could also mimic the dynamics taking place in a complicated phase space, e.g. for glassy systems . From this point of view, sites on the lattice represent microscopic configurations and edges allowed moves from a configuration to another. At low temperatures, most edges correspond to very improbable jumps and may be erased. The tail of density of states of Laplacian on random graphs was studied by means of heuristic arguments by Bray and Rodgers . Localized eigenvectors, closely related to metastable states are of particular relevance for asymptotic dynamics.

Remarkably, the above examples lead to the study of the spectral properties of random symmetric matrices $`𝐖`$ sharing common features. In amorphous media, the elastic energy is a quadratic function of the displacements of the particles from their instantaneous “frozen” positions. The INM are the eigenmodes of the stiffness matrix $`𝐖`$. As for diffusion on random lattices, $`𝐖`$ simply equals the Laplacian operator. In both cases, each row of $`𝐖`$ is comprised of a small (with respect to the size $`N`$ of the matrix) and random number of non-zero coefficients $`W_{ij}`$ and most importantly, diagonal elements fluctuate : $`W_{ii}=_{j(i)}W_{ij}`$ .

In this letter, we present a quantitative approach to explain the spectral properties and especially localization effects of such a random matrix $`𝐖`$ in the simplest case, that is when all off-diagonal elements of $`𝐖`$ are independent random variables. Our analytical approximation is corroborated by exact numerical diagonalizations. We then expose a systematic diagrammatic expansion allowing for the study of more realistic models in presence of correlated $`W_{ij}`$’s.

The spectral properties of $`𝐖`$ can be obtained through the knowlegde of the resolvent $`G(\lambda +iϵ)`$, that is the trace of $`((\lambda +iϵ)\mathrm{𝟏}𝐖)^1`$ . Denoting the average over disorder by $`\overline{()}`$, the mean density of states reads

$$p(\lambda )=\frac{1}{\pi }\underset{ϵ0^+}{lim}\text{Im}\overline{G(\lambda +iϵ)}.$$
(1)   
The averaged resolvent is then written as the propagator of a replicated Gaussian field theory 

$`\overline{G(\lambda +iϵ)}`$ $`=`$ $`\underset{n0}{lim}{\displaystyle \frac{i}{Nn}}{\displaystyle \underset{i}{}d\stackrel{}{\varphi }_i\underset{k=1}{\overset{N}{}}\stackrel{}{\varphi }_k^{\mathrm{\hspace{0.33em}2}}\underset{i}{}z_i\overline{\underset{i<j}{}(1+u_{ij})}}`$ (2)   
$`\text{where}z_i`$ $``$ $`z(\stackrel{}{\varphi }_i)=\mathrm{exp}\left({\displaystyle \frac{i}{2}}(\lambda +iϵ)\stackrel{}{\varphi }_i^{\mathrm{\hspace{0.33em}2}}\right)`$ (3)   
$`u_{ij}`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{i}{2}}W_{ij}(\stackrel{}{\varphi }_i\stackrel{}{\varphi }_j)^2\right)1.`$ (4)   
Replicated fields $`\stackrel{}{\varphi }_i`$ are $`n`$-dimensional vector fields attached to each site $`i`$. To ligthen notations, we have restricted in (4) to the scalar case. We shall focus later on $`𝐖`$ having an internal dimension, as happens to be in the INM problem.

In the uncorrelated case, the $`W_{ij}`$’s ($`i<j`$) are independently drawn from a probability law $`𝒫`$. To take into account geometrical randomness only, we focus on distribution $`𝒫\left(W_{ij}\right)=\left(1\frac{q}{N}\right)\delta \left(W_{ij}\right)+\frac{q}{N}\delta \left(W_{ij}w\right)`$ . Such a bimodal law merely defines a random graph : $`i`$ and $`j`$ can be said to be connected if and only if $`W_{ij}`$ does not vanish. Due to the scaling of the edge-probability $`\frac{q}{N}`$, the mean site-connectivity $`q`$ remains finite for large sizes $`N`$. We rescale the eigenvalues by choosing $`w=\frac{1}{q}`$ to ensure that the support of the spectrum is positive and bounded when $`q\mathrm{}`$ .

Numerical diagonalizations of the random Laplacian $`𝐖`$ have been carried out for different sizes, up to $`N=3200`$. To each eigenvector $`\psi _{i,\mathrm{}}`$ of eigenvalue $`\lambda _{\mathrm{}}`$ normalized to unity is associated the inverse participation ratio $`w_{\mathrm{}}^4=_i|\psi _{i,\mathrm{}}|^4`$. We then define $`w^4(\lambda )d\lambda `$ as the sum of $`w_{\mathrm{}}^4`$ over all $`\psi _{i,\mathrm{}}`$ lying in the range $`\lambda \lambda _{\mathrm{}}\lambda +d\lambda `$, divided by the number $`Np(\lambda )d\lambda `$ of such eigenvectors. Fig. 1 displays $`p(\lambda )`$ and $`w^4(\lambda )`$ for a mean connectivity $`q=20`$ . The central part of the spectrum ($`\lambda _{}<\lambda <\lambda _+`$) has a smooth bell shape and corresponds to extended states. For increasing sizes $`N`$ and at fixed $`\lambda `$, $`w^4`$ vanishes as $`1/N`$ and the breakdown of this scaling identifies the mobility edges : $`\lambda _{}0.47\pm 0.01`$ and $`\lambda _+1.67\pm 0.03`$. Outside the central region, that is for small or large eigenvalues, the eigenstates become localized and the density exhibits successive regular peaks.

We have measured for each eigenvector $`\psi _{i,\mathrm{}}`$ the connectivity $`c_{\mathrm{}}`$ of its center, that is the site $`i_0`$ with maximum component $`|\psi _{i_0,\mathrm{}}|`$. The mean connectivity $`c(\lambda )`$ of the centers of eigenvectors having eigenvalue $`\lambda `$ is plotted fig. 1. It is a smooth monotonous function of $`\lambda `$ in the extended part of the spectrum. In the localized region, $`c(\lambda )`$ is constant over a given peak and integer-valued ($`cc_{}=10`$ on the left side of the spectrum, $`cc_+=33`$ on the right side); the center connectivity abruptly jumps when $`\lambda `$ crosses the borders between peaks. Furthermore, Table 1 shows the good agreement between the weight of peak associated to connectivity $`c`$ and the fraction of sites having $`c`$ neighbors, given by a Poisson law of parameter $`q`$ .

Therefore, numerics indicate that localized eigenvectors are centered on geometrical defects, that is on sites whose number of neighbors is much smaller or much larger than the average connectivity. To support this observation, it is intructive to consider a simpler model including a unique defect i.e. a Cayley tree with connectivity $`c`$ for the central site and $`q+1`$ for all other points  (locally a random graph is equivalent to a tree since no loops of finite length are present). Looking for a localized state $`\psi _i`$ with a radial symmetry $`\psi _i=\psi _{d(i)}`$ where $`d(i)`$ is the distance between site $`i`$ and center, the eigenvalue equations read $`c(\psi _0\psi _1)=q\lambda \psi _0`$ and $`(q+1)\psi _dq\psi _{d+1}\psi _{d1}=q\lambda \psi _d`$ for $`d1`$ . The eigenvalue problem reduces to the search of the solution of a homogeneous linear difference equation of order two (eqn. for $`d1`$) fulfilling a boundary condition (eqn. for $`d=0`$). After a little algebra we have found that strong defects, such that $`|qc|>\sqrt{q}`$ give rise to localized states around the central site with an eigenvalue $`\lambda =\frac{c}{q}(1\frac{1}{qc})`$. The predicted connectivities at mobility edges ($`c_{}=15`$ and $`c_+=25`$ for $`q=20`$) are in poor agreement with numerical findings. A more refined picture requires to take into account the connectivity fluctuations of the neighbours of the central site. We have thus considered a Cayley tree with a coordination number $`c`$ for the central site, $`c^{}`$ for the nearest neighbours and $`q+1`$ for all other points. We have found that localized states due to weak (respectively strong) central connectivity $`c`$ can disappear and become extended if the connectivity of the neighbors $`c^{}`$ reaches large (resp. small) values. In other words, a defect can be screened by an opposite connectivity fluctuation of its surrounding neighbors. Numerics supports this scenario. As $`\lambda `$ varies, $`w^4(\lambda )`$ exhibit oscillations (interpreted as finite-size contributions coming from extended states) of rapidly decreasing amplitudes with increasing $`N`$. These oscillations are correlated (positively for small $`c`$ and negatively for large $`c`$) with the fluctuations of the neighbors connectivity $`c^{}(\lambda )`$ around its mean value $`q+1`$, see inset of fig. 1.

Let us see how the above results may be recovered from theory. Due to the statistical independence of the $`W_{ij}`$’s, the $`u_{ij}`$ (4) interactions are averaged out separately . The resulting theory is of course invariant under any relabelling of the sites $`i`$ and depends on the fields $`\stackrel{}{\varphi }_i`$ through the density $`\rho (\stackrel{}{\varphi })`$ of sites $`i`$ carrying fields $`\stackrel{}{\varphi }_i=\stackrel{}{\varphi }`$ only . The functional order parameter $`\rho (\stackrel{}{\varphi })`$ is found when optimizing the “free-energy” 

$`\mathrm{ln}\mathrm{\Xi }[\rho ]`$ $`=`$ $`{\displaystyle 𝑑\stackrel{}{\varphi }\rho (\stackrel{}{\varphi })\left[\mathrm{ln}z(\stackrel{}{\varphi })\mathrm{ln}\rho (\stackrel{}{\varphi })+1\right]}`$ (5)   
$`+`$ $`{\displaystyle \frac{q}{2}}{\displaystyle 𝑑\stackrel{}{\varphi }𝑑\stackrel{}{\psi }\left(e^{i(\stackrel{}{\varphi }\stackrel{}{\psi })^2/2q}1\right)\rho (\stackrel{}{\varphi })\rho (\stackrel{}{\psi })},`$ (6)   
under the normalization constraint $`𝑑\stackrel{}{\varphi }\rho (\stackrel{}{\varphi })=1`$; $`z(\stackrel{}{\varphi })`$ has been defined in (3). This order parameter is simply related to the original random matrix problem through

$$\rho (\stackrel{}{\varphi })=\frac{1}{N}\overline{\underset{i=1}{\overset{N}{}}C_i\mathrm{exp}\left(\frac{i\stackrel{}{\varphi }^2}{2[(\lambda +iϵ)\mathrm{𝟏}𝐖]_{ii}^1}\right)},$$
(7)   
The $`C_i`$ are normalization constants going to unity as $`n`$ vanishes. Therefore, the averaged resolvent reads $`\overline{G(\lambda +iϵ)}=ilim_{n0}𝑑\stackrel{}{\varphi }\rho (\stackrel{}{\varphi })(\varphi ^1)^2`$.

Finding an exact solution to the maximization equation $`\delta \mathrm{ln}\mathrm{\Xi }/\delta \rho (\stackrel{}{\varphi })=0`$ seems to be a hopeless task. This is a general situation, which arises in the study of the physics of dilute systems (for the case of sparse random matrices see, for example,  ).

Identity (7) may however be used as a starting point for an effective medium approximation (EMA). In the extended part of the spectrum, we expect all matrix elements appearing in (7) to be of the same order of magnitude and thus $`\rho (\stackrel{}{\varphi })`$ to be roughly Gaussian. EMA is therefore implemented by inserting the Gaussian Ansatz

$$\rho ^{EMA}(\stackrel{}{\varphi })=\left(2\pi ig(\lambda )\right)^{\frac{n}{2}}\mathrm{exp}\left(\frac{i\stackrel{}{\varphi }^2}{2g(\lambda )}\right),$$
(8)   
into functional $`\mathrm{\Xi }`$ (6). The average EMA resolvent $`g`$ is then obtained through optimization of $`\mathrm{ln}\mathrm{\Xi }[g(\lambda )]`$. The resulting spectrum, which is given by the imaginary part of $`g`$ divided by minus $`\pi `$, is shown fig. 1. As expected, EMA gives a sensible estimate of the spectral properties in the extended region and of the mobility edges $`\lambda _{}^{EMA}=0.468`$, $`\lambda _+^{EMA}=1.732`$. However, EMA is intrinsically unable to reflect geometry fluctuations and thus the presence of localized states .

To do so, we start by writing the extremization condition of $`\mathrm{ln}\mathrm{\Xi }`$ over $`\rho `$ as

$$\rho (\stackrel{}{\varphi })=[\rho ](\stackrel{}{\varphi }),$$
(9)   
where the functional $``$ may be expanded as

$$[\rho ](\stackrel{}{\varphi })=hz(\stackrel{}{\varphi })\underset{k=0}{\overset{\mathrm{}}{}}\frac{e^qq^k}{k!}\left[𝑑\stackrel{}{\psi }\rho (\stackrel{}{\psi })e^{i(\stackrel{}{\varphi }\stackrel{}{\psi })^2/2q}\right]^k$$
(10)   
$`h`$ is a multiplicative factor equal to unity in the $`n0`$ limit. Equations (9,10) describe an elementary lattice of one central site connected to $`k`$ neighbors according to a Poisson distribution of mean $`q`$. Neighbors carry information about the random matrix elements through the order parameter $`\rho (\stackrel{}{\psi })`$ (7) and interact with the central site via kernel $`\mathrm{exp}(i(\stackrel{}{\varphi }\stackrel{}{\psi })^2/2q)`$ (4,10). Self-consistency requires that the resulting order parameter at central site equals $`\rho `$ . Baring in mind the localization mechanism unveiled in previous paragraphs, we propose a single defect approximation  (SDA). SDA amounts to make the central site interact with $`k`$ neighbours belonging to the effective medium defined above. Since EMA precisely washes out any local geometrical fluctuation, we partially reintroduce them by allowing the connectivity $`k`$ of the central site (the defect) to vary. The SDA order parameter is thus obtained through an iteration of equation (9)

$$\rho ^{SDA}(\stackrel{}{\varphi })=[\rho ^{EMA}](\stackrel{}{\varphi }).$$
(11)   
Using the EMA resolvent $`g`$ (8), we have computed the SDA spectrum for $`q=20`$. The SDA extended part is shown to be in better agreement with numerical results than EMA on fig. 1. Improvement is even more spectacular for localized states that were absent within EMA. We have found Dirac peaks whose weights and eigenvalues are listed Table 1. The agreement with numerical results is quite good. We have verified analytically that SDA peaks do correspond to localized states by calculating $`lim_{ϵ0}[ϵ\overline{G(\lambda +iϵ)G(\lambda iϵ)}]`$  using SDA with two groups of replicas. This quantity gives also access to $`w^4(\lambda )`$, whose value seems slightly higher than numerical measures close to the mobility edges.

Starting from any sensible $`\rho `$, successive iterations of equation (9) would also provide more and more accurate descriptions of the “fractal” structure of the localized peaks at a price of heavier and heavier calculations. Besides being theoretically founded, SDA has the advantage that a single iteration from EMA (which is easily computable) succeeds in capturing localized states in a quantitative way .

In general, the components of $`𝐖`$ are correlated and the average over disorder requires an expansion in terms of the connected correlation functions of the $`u_{ij}`$ interactions

$$\overline{\underset{i<j}{}(1+u_{ij})}=\mathrm{exp}\left(\underset{i<j}{}\overline{u_{ij}}+\frac{1}{2}\underset{i<j,k<l}{\overset{}{}}\overline{u_{ij}u_{kl}}^c+\mathrm{}\right),$$
(12)   
where the prime indicates that the sum runs over different pairs of sites. Free-energy (6) corresponds to the case where all terms in (12) but the first one vanish. The presence of these cumulants (of order $`2,3,\mathrm{}`$ in $`u`$) will result in the addition of cubic, quartic, … $`\rho (\stackrel{}{\varphi })`$ interactions terms to $`\mathrm{ln}\mathrm{\Xi }`$. Though calculations become more difficult, the existence of a variational free-energy $`\mathrm{ln}\mathrm{\Xi }`$ is preserved. This is all what is needed to derive EMA and the optimization equation (11).

Let us see how SDA can be implemented to determine the INM spectrum of amorphous media. We shall restrict to liquids but our approach could also be applied to glasses using the formalism recently developed by Mézard and Parisi . Particles $`i`$ are individuated by their positions $`𝐱_i`$ and interact through a two-body potential $`V(𝐱_i𝐱_j)`$ (hereafter bold letters will denote vectors in the $`D`$-dimensional real space). For a given microscopic configuration, INM are the eigenmodes of the $`D\times N`$-dimensional matrix $`W_{ij}=^2V(𝐱_i𝐱_j)/𝐱_i𝐱_j`$ ($`ij`$). The calculation of the spectrum and the average over particle configurations (with the equilibrium Boltzmann measure at inverse temperature $`\beta `$) can be performed at the same time by introducing a generalized liquid . Each particle is assigned a “position” $`𝐫_i=(𝐱_i,\stackrel{}{\varphi }_i)`$ and the generalized fugacity reads $`z^{}(𝐫_i)=yz(\stackrel{}{\varphi }_i)`$ where $`y`$ is the liquid fugacity and $`z`$ is defined in (3). The grand-canonical partition function $`\mathrm{\Xi }`$ may then diagrammatically expanded in powers of the Mayer bond $`b(𝐫_i,𝐫_j)=\mathrm{exp}(\beta V(𝐱_i𝐱_j)+\frac{i}{2}W_{ij}(\stackrel{}{\varphi }_i\stackrel{}{\varphi }_j)^2)1`$ (the summation over the $`D^2`$ internal indices of $`𝐖`$ is not written explicitely for sake of simplicity). With these notations, $`\mathrm{\Xi }`$ coincides with formula (3.21) of Ref.. It is now straigthforward to take advantage of the variational formulation of the diagrammatic virial expansion by Morita and Hiroike . The generalized density $`\rho (𝐫)=\rho (𝐱,\stackrel{}{\varphi })`$ of particles optimizes

$$\mathrm{ln}\mathrm{\Xi }[\rho ]=𝑑𝐫\rho (𝐫)\left[\mathrm{ln}z^{}(𝐫)\mathrm{ln}\rho (𝐫)+1\right]+𝒮,$$
(13)   
where $`𝒮`$ is the sum of all diagrams composed of bonds $`b(𝐫,𝐫^{})`$ and vertices weighted with $`\rho (𝐫)`$ that cannot be split under the removal of a single vertex, see equation (4.6) of Ref.. Due to translational invariance in real space, $`\rho `$ does not depends on $`𝐱`$ and we are left with a variational functional $`\mathrm{\Xi }`$ of the density $`\rho (\stackrel{}{\varphi })`$. Note that (13) contains (6) as a special case when $`𝒮`$ includes only the simplest single bond diagram. The random graph model we have studied in this letter may be seen as a physical system for which keeping the first coefficient of the virial expansion only is exact.

To our knowledge, Morita and Hiroike’s work has not been used so far in the context of INM theory as a short cut to avoid tedious diagrammatical calculations. In addition, the variational formulation of  allows to implement SDA in a practical way. We are currently attempting to apply the present formalism to characterize localized eigenstates in two- and three-dimensional granular media.

Acknowledgements : We are deeply indebted to D.S. Dean for numerous and thorough discussions on this work. We also thank D.J. Thouless for an enlightening discussion, particularly about the Cayley tree argument.