text:
# Heteroclinic Contours in Neural Ensembles and the Winnerless Competition Principle

## 1 Introduction

Computing with a dynamical system implies that this system changes its behavior depending on the quality and quantity of the incoming information. This is an enormous field and we will concentrate here only on the concept of Winnerless Competition (WLC) that, as we think, is a general principle for information processing by dynamical systems.

Information processing with WLC dynamics is a new area for theoretical study. However, WLC itself has already been observed in many well known experiments in hydrodynamics \[Busse & Heikes, 1980\], population biology \[May & Leonard, 1975\] and laser dynamics \[Roy, 1999\]. For example, the convective roll patterns in a rotating plane layer demonstrate a sequential changing of direction as a result of the competition between different roll orientations (when the rotation rate is large enough, e.g. Kuppers-Lortz instability \[Küppers & Lortz, 1969\]). For a large Prandtle number the critical angle is close to $`60^o`$ (for the steady state rolls it is $`120^o`$). Due to the appearance of new rolls that are also unstable, no pattern becomes a winner and, as a result of the competition, they switch sequentially. Small non-Boussinesq effects are able to make such sequence periodic \[Rabinovich et al., 2000a\].

The discussed process can be described by Lottka-Volterra equations that are well known in population biology and can explain the competition between three spices \[May & Leonard, 1975\]:

$`\dot{a}_1`$ $`=`$ $`a_1[1(a_1+\rho _{12}a_2+\rho _{13}a_3)]`$
$`\dot{a}_2`$ $`=`$ $`a_2[1(a_2+\rho _{21}a_1+\rho _{23}a_3)]`$
$`\dot{a}_3`$ $`=`$ $`a_3[1(a_3+\rho _{31}a_1+\rho _{32}a_3)]`$
In biology the variable $`a_i(t)`$ is the number of individuals in the $`i`$–th population at time $`t`$, and $`\rho _{ij}`$ are competition coefficients measuring the extent to which the $`j`$th species affects the growth rate of the $`i`$–th species. In convection, the variables $`a_i`$ are the intensities of the competitive modes: $`a_i=|c_i|^2`$, and we suppose that the vertical component of the velocity field is (in the limit of small amplitudes) $`u_x=f(z)_{j=1}^3c_j(t)exp\{i𝐤_j𝐫\}`$, where $`z`$ is the component of the position vector $`𝐫`$ in the vertical direction and $`𝐤_j`$ are the wave vectors \[Busse & Heikes, 1980, Rabinovich et al., 2000a\].

The non-symmetry of the coefficients $`\rho _{ij}`$, for example $`\rho _{12}=\rho _{23}=\rho _{31}\rho _+>1`$, $`\rho _{21}=\rho _{32}=\rho _{13}\rho _{}<1`$, guarantees the WLC behavior of the discussed dynamical system. The mathematical image of such behavior is a heteroclinic contour in the phase space $`a_1(t)`$, $`a_2(t)`$, $`a_3(t)`$ (see Fig. 1). Results related to symmetric cases have been reporter earlier \[Ashwin and Field, 1999\] and \[Ashwin and Chossat, 1998\].

The questions that we are going to discuss below are: (i) what are the conditions for the robustness of the WLC, i.e. the topological similarity of the perturbed and original heteroclinic contour; and (ii) how subsystems with WLC behavior interact with each other.

The paper is organized in the following way. First we describe a class of models that use the WLC principle for the representation and processing of incoming information (Sec. 2). Then we discuss the existence and stability of the heteroclinic contour (Sec. 3), and the robustness, i.e. birth of a stable limit cycle in the vicinity of the destroyed heteroclinic loop in a perturbed system (Sec. 4) with some examples from computer modeling. Finally we discuss several WLC strategies used by living neural systems to perform complex information processing (Sec. 5).

## 2 The Models

The activity of many different neural networks \[Rabinovich et al., 2001, Varona et al., 2001, Abeles et al., 1995, Cohen & Grossberg, 1983\], can be described qualitatively with the following dynamics:

$`\dot{a}_i`$ $`=`$ $`a_i(\sigma (𝑯,𝑺){\displaystyle \underset{j=1}{\overset{N}{}}}\rho _{ij}a_j+H_i(t))+S_i(t)`$ (1)
where $`a_i>0`$ represents the instantaneous spiking rate of the principal neurons (PNs) that are making the computation, $`\rho _{ij}`$, represents the strength of inhibition in $`i`$ by $`j`$, $`H_i(t)`$ represents the action from other neural ensembles, and $`S_i(t)`$ represents the stimuli from the sensors. In many neural networks, the inhibition among PNs is the result of the action of inhibitory local neurons (LNs). Usually LNs also receive an external input and because of this $`\rho _{ij}`$ can depend on the stimuli.

The dynamical system (1) in the case $`\sigma =1`$, $`𝑯(t)=𝑺(t)=0`$ is the Lottka-Volterra model. The dynamics of the system is well known when the matrix $`\rho _{ij}`$ is symmetric ($`\rho _{ij}=\rho _{ji}`$). In this case the autonomous system has a global Lyapunov function \[Cohen & Grossberg, 1983, Hopfield, 1982\] and every trajectory approaches one of the numerous possible equilibrium points. For example, if the inhibitory connections are identical, $`\rho _{ij}=\rho `$, $`\rho _{ii}=1`$, this system has only one global attractor, e.g. $`a_i=a_0=1/[1+\rho (N1)]`$ for $`\rho <1`$, and $`N`$ attractors: $`a_i=a_0=1`$, $`a_{ji}=0`$ if $`\rho >1`$. No other attractors, e.g. limit cycles, or strange attractors are present in the system. The situation is much more complex and interesting when the inhibition is non-symmetric: $`\rho _{ij}\rho _{ji}`$. A detailed analysis is only possible in the case $`N=3`$ (see references \[Afraimovich et al., 2001, Chi et al., 1998, Rabinovich et al., 2001\]). When $`\rho _{ij}>1,\rho _{ji}<1`$ there exists a heteroclinic contour in the phase space of the system that consists of saddle points and one-dimensional separatrices connecting them. In some regions of the parameter space $`\{\rho _{ij}\}`$, such heteroclinic contour (or limit cycle in its vicinity) is a global attractor. If $`\rho _{ij}`$ depend on the stimulus, e.g. as a result of a learning mechanism, the system (1) can generate different heteroclinic contours for different stimuli \[Rabinovich et al., 2001\].

Suppose the matrix

$`(\rho _{ij})=\left(\begin{array}{ccc}1& \alpha _1& \beta _1\\ \beta _2& 1& \alpha _2\\ \alpha _3& \beta _3& 1\end{array}\right)`$

and $`0<\alpha _i<1<\beta _i`$ and $`\kappa _i=(\beta _i1)/(1\alpha _i)`$. Then the heteroclinic contour is a global attractor if $`\kappa _1\kappa _2\kappa _3>1`$, and the nontrivial fixed point $`A(a_1^0,a_2^0,a_3^0)`$ is a saddle point. If $`\kappa _1\kappa _2\kappa _3=1`$, this fixed point becomes neutrally stable and there exists a family of neutrally stable periodic solutions in the phase space. When $`\kappa _1\kappa _2\kappa _3<1`$, $`A`$ becomes a global attractor. The heteroclinic orbit exists but looses its stability. It is important to emphasize that in the case $`\kappa _1\kappa _2\kappa _3>1`$ a small perturbation is able to destroy the heteroclinic orbit and then a stable limit cycle appears in its vicinity. This limit cycle is characterized by a finite time period of switching among different states, in contrast with the infinite time of motion along the heteroclinic loop.

When $`N>3`$ the dynamics of system (1) can be very complex and even chaotic \[Varona et al., 2001\]. Here we are interested in the existence and stability of the heteroclinic contours, which are the mathematical image of the winnerless competition behavior. Such orbits may exist only in the nonsymmetric case e.g. $`\rho _{ij}\rho _{ji}`$, when saddle points (in the heteroclinic contours) satisfy several conditions.

## 3 Existence and Stability of the Heteroclinic Contour

In this section we consider the canonical Lottka-Volterra model

$`\dot{a}_i`$ $`=`$ $`a_i[1(a_i+{\displaystyle \underset{ij}{\overset{N}{}}}\rho _{ij}a_j)],`$ (2)
and derive conditions of existence and stability of heteroclinic contours.

### 3.1 Necessary Conditions

#### 3.1.1 ”Codimension one” saddle points

A heteroclinic contour consists of finitely many saddle equilibria and finitely many heteroclinic orbits connecting these equilibria. Let’s denote by $`A_1`$ the equilibrium point $`(1,0,0,\mathrm{},0)`$, by $`A_2`$ the point $`(0,1,\mathrm{}0)`$, and by $`A_N`$ the point $`(0,0,\mathrm{},1)`$. For the sake of simplicity we assume that there is a heteroclinic orbit $`\mathrm{\Gamma }_{ii+1}`$ connecting the points $`A_i`$ and $`A_{i+1}`$, $`i=1,\mathrm{}N`$ and $`A_{N+1}A_1`$. (If not, we can always apply a change of variables in the form of a permutation). The contour can serve as an attracting set if every point $`A_i`$ has only one unstable direction. By direct verification it can be shown that $`A_i`$ satisfies this assumption provided that:

$`\rho _{ki}>1,ki+1,`$ (3)
and

$`\rho _{i+1i}<1.`$ (4)
(Here $`i+1=1`$ if $`i=N`$).

Moreover, if (3) and (4) are satisfied then the unstable direction at the point $`A_i`$ is parallel (at that point) to the ort $`(0\mathrm{}010\mathrm{}0)`$, where 1 corresponds to the $`i`$–th coordinate. An intersection of hyperplanes $`P_{2i}=_{j=1,ji,i+1}^N\{a_j=0\}`$ is a two-dimensional invariant manifold containing points $`A_i`$ and $`A_{i+1}`$ such that $`A_i`$ is a saddle point on $`P_{2i}`$ and $`A_{i+1}`$ is a stable node on $`P_{2i}`$. The system (2) on $`P_{2i}`$ has the form:

$`\dot{a}_i`$ $`=`$ $`a_i[1(a_i+\rho _{ii+1}a_{i+1})]`$ (5)
$`\dot{a}_{i+1}`$ $`=`$ $`a_{i+1}[1(a_{i+1}+\rho _{i+1i}a_i)]`$ (6)
and, from (3) and (4), one has $`\rho _{ii+1}>1`$, $`\rho _{i+1i}<1`$.

This implies that there are no equilibrium points in the region $`a_i>0`$, $`a_{i+1}>0`$, and since $`\dot{a}_{i+1}<0`$ if $`a_{i+1}1`$ then it is simple to see that the separatrix, say $`\mathrm{\Gamma }_i`$ of the saddle point $`A_i`$ must go to the attractor $`A_{i+1}`$, i.e. there is a heteroclinic connection between $`A_i`$ and $`A_{i+1}`$ on the plane $`P_{2i}`$ (For the case N=3 see the proof in \[Waltman, 1983\]).

#### 3.1.2 Leading directions

The point $`A_i`$ on $`P_{2i}`$ is a stable node with characteristic numbers $`\lambda _1=1`$ and $`\lambda _2=1\rho _{ii+1}`$. The leading direction at $`A_{i+1}`$ is determined by the absolute values of $`\lambda _1`$ and $`\lambda _2`$: if $`\lambda _1>\lambda _2`$ then the leading direction is parallel to the $`a_{i+1}`$-axis, and if $`\lambda _1<\lambda _2`$ then the leading direction is transversal to the $`a_{i+1}`$-axis on $`P_{2i}`$. We assume that the last inequality holds, i.e.

$$\rho _{ii+1}<2,$$
(7)
then the majority of orbits, (including $`\mathrm{\Gamma }_i`$) go to $`A_{i+1}`$ following a direction $`\stackrel{}{l}=(1,\rho _{i+1i}/(2\rho _{ii+1}))`$ transversal to the $`a_{i+1}`$-axis.

The vector $`\stackrel{}{l}`$ on $`P_{2i}`$ can be embedded into the hyperplane $`H_i:\{a_{i+2}=0\}`$ as $`\stackrel{}{L}=(00\mathrm{}1,\rho _{i+1i}/(2\rho _{ii+1}),0\mathrm{}0)`$, $`1`$ on the $`i`$–th place, and one can ask if the direction $`\stackrel{}{L}`$ is the leading direction for the node $`A_{i+1}`$ on this hyperplane. To see sufficient conditions for that assumption we have to take into account that the characteristic numbers at point $`A_{i+1}`$ of the system (2) restricted to the hyperplane $`\{a_{i+2}=0\}`$ are $`1\rho _{ii+1},\mathrm{},1\rho _{i1i+1},1,1\rho _{i+2i+1},\mathrm{},1\rho _{Ni+1}`$ (they are all negative because of (3)). Hence, if

$`\rho _{ki+1}>\rho _{ii+1},ki,`$ (8)
then $`1\rho _{ii+1}`$ is the characteristic value closest to zero and $`\stackrel{}{L}`$ is the leading direction at $`A_{i+1}`$ on $`H_i`$. We assume that (8) is satisfied. This condition is not necessary for the validity of the results below, but it essentially simplifies the description of the results and calculations.

#### 3.1.3 Saddle values

The point $`A_i`$ is a saddle point on $`P_{2i}`$. One can write a map from a transversal to the stable separatrix into a transversal to the unstable separatrix along the orbits going through a neighborhood of $`A_i`$ (see for instance \[Shilnikov et al., 2001\]). In suitable coordinates $`(\xi ,\eta )`$ it has the form:

$`\xi =c\eta ^{\nu _i}`$ (9)
where $`\eta `$ is a deviation from the stable manifold, $`\xi `$ is a deviation from the unstable one, $`c`$ is a constant and

$`\nu _i={\displaystyle \frac{1\rho _{ii+1}}{1\rho _{i+1i}}}{\displaystyle \frac{\rho _{ii+1}1}{1\rho _{i+1i}}}`$ (10)
is the “saddle value” (\[Shilnikov et al., 2001\]). If $`\nu _i>1`$ then the map (9) is a local contraction and $`P_i`$ is a dissipative saddle. If $`\nu _i<1`$ then (9) is a local expansion.

#### 3.1.4 Stability of the heteroclinic contour

The following result tells us that the contour $`\mathrm{\Gamma }=_{i=1}^N\mathrm{\Gamma }_iA_i`$ can be an attractor.

###### Theorem 1

Assume that conditions (3), (4), (7), (8) are satisfied and

$`\nu ={\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{\rho _{ii+1}1}{1\rho _{i+1i}}}>1`$ (11)
(here $`i+1=1`$ if $`i=N`$). Then there is a neighborhood $`U`$ of the contour $`\mathrm{\Gamma }`$ such that for any initial condition $`a^0=(a_1^0,\mathrm{},a_N^0)`$ in $`U`$ with $`a_i^0>0`$, one has dist$`(a(t),\mathrm{\Gamma })0`$ as $`t\mathrm{}`$ where $`a(t)`$ is the orbit going through $`a^0`$.

### 3.2 Proof of Theorem 1

The proof of the theorem is based on the construction of the Poincaré map along orbits in a neighborhood of the contour $`\mathrm{\Gamma }`$. Let $`W_i^s`$ ($`W_i^u`$) be a stable (unstable) manifold of the point $`A_i`$ and $`P_i`$ ($`Q_i`$) be a point on the heteroclinic orbit $`\mathrm{\Gamma }_{i1}`$ ($`\mathrm{\Gamma }_i`$) in a small neighborhood of $`A_i`$–see Fig 2.

#### 3.2.1 Local map

Let $`S_{P_i}`$ ($`S_{Q_i}`$) be a piece of a transversal to $`\mathrm{\Gamma }_{i1}`$ (to $`\mathrm{\Gamma }_i`$) hyperplane going through $`P_i`$ (through $`Q_i`$). Without loss of generality, we may assume that $`S_{P_i}`$ is a piece of a hyperplane, parallel to the hyperplanes $`a_{i+1=0}`$. A local map $`f_i:S_{P_i}S_{Q_i}`$ along orbits in a neighborhood of $`A_i`$ is well defined. In suitable coordinates it has the form \[Shilnikov et al., 2001\].

$`\xi _i=c_i\eta _i^{\nu _i},\chi _i=\phi _i(y_i,\eta _i)`$ (12)
where $`\eta _i\mathrm{}`$ is a coordinate on $`S_{P_i}`$, “parallel” to $`W_i^u`$, $`y\mathrm{}^{N2}`$ is a vector-coordinate transversal to the $`\eta _i`$-axis on $`S_{P_i}`$, $`\xi _i`$ is a coordinate on $`S_{Q_i}`$, “parallel” to the leading direction on $`W_i^S`$ at $`A_i`$, $`y_i`$ is a vector-coordinate transversal to the $`\xi _i`$-axis on $`S_{Q_i}`$. Moreover,

$`|{\displaystyle \frac{\phi _i}{y_i}}|\overline{c}_i\eta _i^{\beta _i},|{\displaystyle \frac{\phi _i}{\eta _i}}|\overline{c}_i|y_i||\eta _i|^{\beta _i1}`$ (13)
where $`\overline{c}_i>0`$ is a constant and

$`\beta _i>\nu _i`$ (14)
#### 3.2.2 Global map

The heteroclinic orbit $`\mathrm{\Gamma }_i`$ has a piece joining the points $`Q_i`$ and $`P_{i+1}`$. Thus, the global map $`F_i:\stackrel{~}{S}_{Q_i}S_{P_{i+1}}`$ along orbits in a neighborhood of this piece is well-defined where $`\stackrel{~}{S}_{Q_i}S_{Q_i}`$ is a small neighborhood of the point $`Q_i`$ on $`S_{Q_i}`$. This map is a diffeomorphism and has the form

$`\eta _{i+1}`$ $`=`$ $`a_{i1}\xi _i+a_{i2}\xi _i+\mathrm{}`$
$`y_{i+1}`$ $`=`$ $`b_{i0}+b_{i1}\xi _i+b_{i2}\xi _i+\mathrm{},`$ (15)
where the dots denote nonlinear terms. The orbit $`\mathrm{\Gamma }_i`$ belongs to the intersection of invariant hyperplanes $`\{a_j=0\}`$, $`ji`$, $`ji+1`$. Therefore the hyperplane $`\{\xi _i=0\}`$ on $`S_{Q_i}`$ is mapped by $`F_i`$ into the hyperplane $`\{\eta _{i+1}=0\}`$ on $`S_{P_{i+1}}`$ which means that $`a_{i2}=0`$ in (15), and

$`a_{i1}0`$ (16)
Hence (15) has the form:

$`\eta _{i+1}`$ $`=`$ $`a_i\xi _i+\mathrm{}`$
$`y_{i+1}`$ $`=`$ $`b_{i0}+b_{i1}\xi _i+b_{i2}\xi _i+\mathrm{},`$ (17)
$`a_{i1}=:a_i`$.

#### 3.2.3 Poincaré map

We may construct a Poincaré map $`F=S_{P_1}S_{P_1}`$ as the superposition of maps $`f_i`$, $`F_i`$, i.e. $`F`$ has the form $`F=F_Nf_N\mathrm{}F_2f_2F_1f_1`$.

The map $`F_if_i`$ has the form

$`\eta _{i+1}`$ $`=`$ $`a_ic_i\eta _i^{\nu _i}+\mathrm{}`$
$`y_{i+1}`$ $`=`$ $`b_{i0}+b_{i1}c_i\eta _i^{\nu _i}+b_{i2}\phi _i(y_i,\eta _i)+\mathrm{}`$ (18)
So, in this approximation the form of the map $`F`$ along the $`\eta `$-coordinates is independent of the $`y`$-coordinates , and one may consider the one-dimensional approximation

$`\eta _{i+1}`$ $`=`$ $`a_ic_i\eta _i^{\nu _i},i=1,\mathrm{},N`$ (19)
$`\eta _{N+1}`$ $`:=`$ $`\overline{\eta }_1\stackrel{~}{F}(\eta _1)`$
It is simple to see that the map $`\stackrel{~}{F}:\eta _1\overline{\eta }_1`$ is a contraction provided that condition (11) is satisfied. Indeed, it follows from (19) that $`\overline{\eta _1}=A\eta _1^\nu `$, where $`A`$ is a constant and $`\nu =_{i=1}^N\nu _i`$. Hence if $`\nu >1`$, then $`\overline{\eta _1}/\eta _1<1`$ if $`\eta _1`$ is small enough, the map (19) is a contraction and $`\eta _1=0`$ is an attracting fixed point (corresponding to the contour $`\mathrm{\Gamma }`$).

Consider now the map $`F_if_i`$. The differential

$`DF_if_i=\left(\begin{array}{cc}a_ic_i\nu _i\eta _i^{\nu _i1}+\mathrm{}\hfill & 0+\mathrm{}\hfill \\ b_{i1}c_i\nu _i\eta _i^{\nu _i1}+b_{i2}\frac{\phi _i}{\eta _i}+\mathrm{}\hfill & \frac{\phi _i}{y_i}+\mathrm{}\hfill \end{array}\right)`$ (22)
and because of (13), (14),

$`DF_if_iB_i\eta _i^{\nu _i1}`$ (23)
where $`B_i`$ is a constant. Since $`DF=_{i=1}^ND(F_if_i)`$, then

$`DFB{\displaystyle \underset{i=1}{\overset{N}{}}}\eta _i^{\nu _i1}`$ (24)
where $`B=_{i=1}^NB_i`$. We estimate now $`_{i=1}^N\eta _i^{\nu _i1}`$ by using (18). In the expression $`\eta _N^{\nu _N1}\eta _{N1}^{\nu _{N1}1}\mathrm{}\eta _2^{\nu _21}\eta _1^{\nu _11}`$ let us estimate first $`\eta _2^{\nu _21}\eta _1^{\nu _11}`$. Since $`\eta _2=a_1c_1\eta _1^{\nu _1}+\mathrm{}`$, then $`\eta _2^{\nu _21}\eta _1^{\nu _11}=(a_1c_1)^{\nu _21}(\eta _1^{\nu _1}+\mathrm{})^{\nu _21}\eta _1^{\nu _11}\text{c}onst\eta _1^{\nu _1\nu _21}`$.

Now $`\eta _3^{\nu _31}\eta _2^{\nu _21}\eta _1^{\nu _11}=(a_2c_2)^{\nu _31}(\eta _2^{\nu _2}+\mathrm{})^{\nu _31}\eta _1^{\nu _11}const\eta _2^{\nu _3\nu _21}\eta _1^{\nu _11}=\text{c}onst(a_1c_1)^{\nu _3\nu _21}(\eta _1^{\nu _1}+\mathrm{})^{\nu _3\nu _21}\eta _1^{\nu _11}\text{c}onst\eta _1^{\nu _3\nu _2\nu _11}`$.

Assume (inductively) that $`\eta _k^{\nu _k1}\mathrm{}\eta _2^{\nu _21}\text{c}onst\eta _2^{\nu _k\mathrm{}\nu _21}`$, then $`\eta _k^{\nu _k1}\mathrm{}\eta _2^{\nu _21}\eta _1^{\nu _11}\text{c}onst(a_1c_1\eta _1^{\nu _1}+\mathrm{})^{\nu _k\mathrm{}\nu _21}\eta _1^{\nu _11}\text{c}onst\eta ^{_{i=1}^k\nu _i1}`$

and therefore,

$`{\displaystyle \underset{i=1}{\overset{N}{}}}\eta _i^{\nu _i1}C\eta _1^{\nu 1}`$ (25)
where $`C`$ is a constant. Hence,

$`DFBC\eta _1^{\nu 1},`$ (26)
i.e. $`F`$ is a contraction in a neighborhood of the point $`\eta _1=0`$, $`y_1=0`$ corresponding to the contour $`\mathrm{\Gamma }`$.

This concludes the proof of Theorem 1.

### 3.3 Selecting the Contour

In this subsection we show how the canonical system (2) can be obtained from a general system for a neural network with inhibitory connections. Suppose that the dynamics of a network of $`N+M`$ inhibitory neurons with dynamical variables $`y_i(t)=(y_i^{(1)}(t),\mathrm{},y_i^{(m)}(t)),i=1,\mathrm{},N+M`$ can be described in the form of the following system of ODEs:

$`\dot{y}_i`$ $`=`$ $`F(y_i){\displaystyle \underset{j=1}{\overset{N+M}{}}}G_{ij}(s)(y_i,y_j)+\stackrel{~}{S}_i(t)`$ (27)
where $`F`$ is a nonlinear function that describes dynamics of individual neuron, $`G_{ij}(S)(y_i,y_j)`$ is a nonlinear operator describing an inhibitory action of the $`j`$-th neuron onto the $`i`$–th neuron, $`S(t)=(S_1(t),\mathrm{},S_{N+M}(t))`$ and $`(\stackrel{~}{S}_1(t),\mathrm{},\stackrel{~}{S}_{N+M}(t))`$ are the vectors representing stimuli to the network. We restrict ourselves to a subclass of systems (27) represented in the form (1). It is known (see above) that a stimulus acts in two ways: (i) it adds the perturbation $`\stackrel{~}{S}(t)`$ into (27) as an external force, and (ii) it forms the matrix $`G_{ij}(S)`$. A simplified model that describes the firing rate of the neurons can be written in the form (1\[Rabinovich et al., 2001\], where $`\sigma =1`$ when there is no stimulus, and $`\sigma =1`$ when the stimulus has a component at neuron $`i`$. In the absence of the external force $`\stackrel{~}{S}(t)`$ the system (2) is just a subsystem of (1) for which all $`\sigma =+1`$.

We describe now an algorithm to obtain the system (2) from (1) (provided that $`\stackrel{~}{S}(t)=0`$). First of all we single out only indexes $`i`$ for which $`\sigma _i(S)=+1`$. As a result we obtain a system which has the form:

$`\dot{a}_i`$ $`=`$ $`a_i[1(a_i+{\displaystyle \underset{ij}{\overset{N+M_1}{}}}\rho _{ij}a_j)]`$ (28)
with $`M_1M`$. Then we consider the only those indexes $`i`$ for which conditions (3) and (4) are satisfied. After some permutation we have the system

$`\dot{a}_i`$ $`=`$ $`a_i[1(a_i+{\displaystyle \underset{ij}{\overset{N+M_2}{}}}\rho _{ij}a_j)]`$ (29)
with $`M_2M_1`$. Let us form a graph $`G`$ now: its nodes are equilibrium points of (29), e.g. points $`A_i=(0\mathrm{}0i=10\mathrm{})`$. There is an edge starting at the point $`A_i`$ and ending at the point $`A_j`$ if there is a separatrix (a piece of $`W_i`$) joining and $`A_i`$ and $`A_j`$. Thus, for any point $`A_i`$ there is the only one edge starting at $`A_j`$. It implies that there exists a subgraph of this graph in the form of a cycle that contains, say, $`N`$ vertices. The equation describing dynamics of $`N`$ corresponding neurons has exactly the form (2) (maybe after some permutative change of variables).

## 4 Birth of a Stable Limit Cycle

A direct corollary of Theorem 1 is the possibility of the birth of the stable limit cycle in system (2) when perturbed in an appropriate way. A perturbation should act in such a way that a Poincaré map $`S_{P_1}S_{P_1}`$ for the perturbed system has an absorbing region and a fixed point inside it. Such a condition can be expressed differently. Let us do it as follows. Consider the system

$`\dot{a}_i=a_i[1(a_i+{\displaystyle \underset{ji}{\overset{N}{}}}\rho _{ij}a_j)]+ϵ\mathrm{\Psi }_i(a)`$ (30)
that coincides with (2) for $`ϵ=0`$, where $`a=(a_1,\mathrm{},a_N)`$ and $`\mathrm{\Psi }_i`$ is a smooth function, $`i=1,\mathrm{},N`$. For small $`ϵ>0`$ the system (30) has saddle equilibrium points $`A_{iϵ}`$ and separatrices $`\mathrm{\Gamma }_{iϵ}`$ (the half of $`W_{iϵ}^u`$ such that $`A_{iϵ}A_i`$, as $`ϵ0`$ and $`lt_{ϵ0}\mathrm{\Gamma }_{iϵ}\mathrm{\Gamma }_{ii+1}`$ (here $`lt`$ means the topological limit, i.e. the set of the accumulation points).

###### Theorem 2

Assume that the conditions of Theorem 1 are satisfied,

$`lt_{ϵ0}({\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{\Gamma }_{iϵ})=\mathrm{\Gamma }`$ (31)
and at least one of the separatrices $`\mathrm{\Gamma }_{iϵ}`$ is not a heteroclinic orbit. Then for any sufficiently small $`ϵ>0`$ the system (30) has a stable limit cycle $`L_ϵ`$ (in a neighborhood of $`\mathrm{\Gamma }`$) such that $`lt_{ϵ0}L_ϵ=\mathrm{\Gamma }`$.

The proof of this Theorem can be done in the standard way, i.e., by construction of the Poincaré map and by showing (as in the proof of Theorem 1) that this map is a contraction in an absorbing region. The condition (31) (or a similar condition) is necessary and sufficient for the existence of an absorbing region. We omit the proof in the present work, since it is really a standard one (the corresponding technique can be found in \[Shilnikov et al., 2001, Afraimovich et al., 1998, Afraimovich et al., 2001\]. Thus, one may say that the system under consideration is robust in the following sense: the attractor of a perturbed system remains in a small neighborhood of the “unperturbed” attractor.

Numerical results show that the system (30) where $`\mathrm{\Psi }_i(a)0`$ satisfies the condition (31) and has a stable limit cycle. An example is shown in Fig. 3. In particular, it is clear from this figure that the system has the only one -global- attractor. In this example, the simulations were performed with the following equations:

$`\dot{a_i}`$ $`=`$ $`a_i(1{\displaystyle \underset{j=1}{\overset{N=6}{}}}\rho _{ij}a_j)+ϵa_ia_{i+3}`$ (32)
where $`i=1,2,\mathrm{},6`$ and $`i+3i3`$ if $`i>3`$. We used the following values of the connection matrix $`\rho _{ij}0`$:

$`\rho _{1,3}`$ $`=`$ $`\rho _{3,5}=\rho _{5,1}=5;\rho _{4,6}=\rho _{2,4}=\rho _{6,2}=2`$
$`\rho _{1,6}`$ $`=`$ $`\rho _{2,1}=\rho _{3,2}=\rho _{4,3}=\rho _{5,4}=\rho _{6,5}=1.5`$
$`\rho _{1,1}`$ $`=`$ $`\rho _{2,2}=\rho _{3,3}=\rho _{4,4}=\rho _{5,5}=\rho _{6,6}=1`$ (33)
with $`ϵ=0.01`$.

Time series showing the switching of activities displayed in this system by the WLC are depicted in Fig. 4.

## 5 WLC in Real Neural Systems

In this section we describe two examples of living neural systems that appear to use WLC strategies to process sensory information. As we discussed in the Introduction, sensory systems need both robustness and sensitivity to acquire and process incoming signals. The WLC among sensory neurons guarantees the coexistence of these essential characteristics and can explain the experimental recordings in several sensory systems. We will illustrate two of these sensory systems using models with WLC dynamics.

### 5.1 WLC in Olfactory Processing

Observed features of olfactory processing networks \[Wehr & Laurent, 1996, Laurent et al., 2001\] can guide the study of computation using competitive networks. In Fig. 5 we show the simultaneously recorded activity of three different projection neurons (PNs) in the locust olfactory system, i.e. antennal lobe (AL), evoked by two different odors: despite similar PN activities before the stimulus onset (the result of the action of noise) each odor evokes a specific spatio-temporal activity pattern that results from interactions between these and other neurons in the network \[Laurent et al., 2001, Bazhenov et al., 2001, Rabinovich et al., 2000b\].

As we have discussed, WLC networks produce identity-temporal or spatio-temporal coding in the form of deterministic trajectories moving along heteroclinic orbits that connect saddle fixed points or saddle limit cycles (see Fig. 1) in the system’s state space. These saddle states correspond to the activity of specific neurons or groups of neurons and the separatrices connecting these states correspond to sequential switching from one state to another.

From the experimental results \[Wehr & Laurent, 1996, Laurent et al., 2001\] we infer that a stimulus acts in two principal ways: (1) it excites a subset of projector neurons; (2) it modifies the effective inhibitory connections between the projector neurons as a result of activation of the inhibitory interneurons that connect different PNs. The intrinsic dynamics of these neurons is governed by many variables corresponding to ion channels and intracellular processes. Such detailed description however is not needed to illustrate the principle of “coding with separatrices”. We need only to capture the ‘firing’ or ‘not-firing’ state of the component neurons. We thus simplify our model to an equation for the firing rate $`a_i(t)>0`$ of neural activity \[Rabinovich et al., 2000b\]:

$$\dot{a}_i(t)=a_i(t)\left[\sigma _i(𝐒)\left(a_i+\underset{ji}{\overset{N}{}}\rho _{ij}(𝐒)a_j(t)\right)\right]+S_i(t).$$
(34)
where $`\rho _{ij}(𝐒)`$ is the strength of inhibition of neuron $`j`$ onto $`i`$. $`\sigma _i(𝐒)=1`$, when there is no stimulus, and $`\sigma (𝐒)=+1`$ when the stimulus has a component at neuron $`i`$. When $`\sigma _i=1`$, the quiet resting state $`a_i=0`$ is stable. When a stimulus is applied and $`\sigma _i=+1`$, the system moves away from this quiet state onto a sequence of heteroclinic trajectories. This instability triggers the system into rapid action, provides robustness against noise and allows a response independent of the state at stimulus onset.

When the inhibitory connections are not symmetric, the system with $`N`$ competitive neurons has different closed heteroclinic orbits that consist of saddle points and one dimensional separatrices connecting them. Such heteroclinic orbits are global attractors in phase space and are found in various regimes of the $`\rho _{ij}(𝐒)`$. This implies that if the stimulus is changed, another orbit in the vicinity of the heteroclinic orbit becomes a global attractor for this stimulus. Such rich behavior can be illustrated also by an inhibitory ensemble of spiking neurons. We have studied a network of inhibitory connected FitzHugh-Nagumo neurons ($`i=1,2,\mathrm{},9`$\[Rabinovich et al., 2000b\]:

$`\tau _1{\displaystyle \frac{dx_i(t)}{dt}}`$ $`=`$ $`f(x_i(t))y_i(t)z_i(t)(x_i(t)\nu )+0.35+S_i`$
$`{\displaystyle \frac{dy_i(t)}{dt}}`$ $`=`$ $`x_i(t)by_i(t)+a`$
$`\tau _2{\displaystyle \frac{dz_i(t)}{dt}}`$ $`=`$ $`{\displaystyle \underset{j}{}}g_{ji}G(x_j(t))z_i(t)`$ (35)
Here we use a dynamical model of inhibition: $`z_i(t)`$ is a synaptic current modeled by first order kinetics. The variable $`x_i(t)`$ denotes the membrane potential, $`y_i(t)`$ is a recovery variable, and $`f(x)=x\frac{1}{3}x^3`$ is the internal FN nonlinearity. The stimulus is taken as a constant. We use a step function for $`G(x)=0,x0,\text{and}G(x)=1,x>0`$, as the synaptic connection. $`S_i`$ is the stimulus, and $`g_{ji}`$, the strength of synaptic inhibition: $`g_{ji}=2`$ if the $`j^{th}`$ neuron inhibits the $`i^{th}`$; $`0`$ otherwise. The other parameters are $`a=0.7,b=0.8,\tau _1=0.08,\tau _2=3.1,\nu =1.5`$.

Our numerical simulations show that the network produces different spatio-temporal patterns in response to different stimuli. Fig. 6 presents examples of these activities corresponding to two different stimuli. The system was in the resting state $`x_i1.2,y_i0.62,z_i=0`$ before the stimulus began at $`t=0`$. As one can see, the patterns are considerably different and distinguishable. The heteroclinic contour in this network consists of a finite number of saddle limit cycles and the same number of heteroclinic orbits connecting these cycles (see, as an example, Fig. 1, right panel). A detailed characterization of this network as an information processing device has been reported in \[Rabinovich et al., 2001\].

### 5.2 WLC in the Gravimetric Neurons of the Mollusk Clione

Neural networks with WLC dynamics are able to generate new information to answer a simple external signal. Such information can be used for the organization of complex activity and, in particular, chaotic behavior of some animals. Let us consider now the hunting activity of a marine mollusk Clione. This mollusk is a predator lacking a visual system. It feeds on a small mollusk, Limacina. The hunting behavior is a random search for prey: Clione ”scans” the surrounding space in order to locate and catch the prey. Such behavior is turned on by the smell of the Limacina. The main role in the organization of such motion of Clione is played by a sensory neural network inside the gravimetric organs: the statocysts (see Fig.7). These special sensory organs are responsible for the orientation in the gravitational field \[Arshavsky et al., 1993\].

It is well known from the physiological data that the statocysts have up to 12 receptor neurons (SRNs) that are coupled with inhibitory synapses \[Arshavsky et al., 1993\]. These neurons respond to the pressure exerted by the statolith, a stone located inside the statocyst. If no information about a prey (received by the chemical receptors) is present, the receptor neuron $`D`$ (down, see Fig. 7) is excited by the statolith and it inhibits other SRNs, and the network responds in a winner-take-all mode. As a result, the information generated by $`D`$ SRN arrives to the corresponding Central Pattern Generators (CPGs) that control the tail and wing movements. These CPGs establish the habitual ”head up” position of Clione’s body. However if a special Hunting Central Neuron (HCN) receives a message from the chemo-sensors about the presence of a prey, HCN excites the SRNs organizing a WLC among them as we will illustrate with a model. The behavior of the Clione in this case does not depend on the direction of the gravitational field and it moves in a random-like trajectory.

For the phenomenological modeling of the statocyst “hunting” dynamics we can neglect the statolith inertial dynamics and take into account the only key point: the position of the mollusk’s body uniquely depends on the message that SRNs are sending to the central neurons that produce the commands to the CPGs. Thus, as a starting point, we consider just a SRN network under the action of the HCN excitation. We suppose that, as a result of the HCN stimulation, all SRNs (”left”, ”right”, ”back”, ”front”, ”down”, and ”up”) are in the same situation: they receive and send two inhibitory synapses (see Fig. 7, right panel).

The dynamics of the SRN’s network can be described by model (1) with $`N=6`$. In this case, $`a_i>0`$ represents the instantaneous spiking rate of the receptor neuron $`i`$, $`H_i(t)`$ represents the stimulus from the hunting neuron to neuron $`i`$, and $`S_i(t)`$ represents the action of the statolith on the receptor that is pressing. When there is no stimulus from the hunting neuron ($`H_i=0,i`$) or the statolith ($`S_i=0,i`$), then $`\sigma (𝑯,𝑺)=1`$ and all neurons are silent; $`\sigma (𝑯,𝑺)=1`$ when the hunting neuron is active and/or the statolith is pressing one of the receptors. In our simulations, we have used the values $`\rho _{ij}0`$ specified in (33).

When there is no activation of the sensory neurons from the hunting neuron, the effect of the statolith ($`S_i0`$) in this model is to induce a higher rate of activity on one of the neurons (the neuron $`i`$ where it rests for a big enough $`S_i`$ value). We assume that this higher rate of activity affects the behavior of the motoneurons to organize the head up position. The other neurons are either silent or have a lower rate of activity and we can suppose that they do not influence the posture of Clione.

When the hunting neuron is active a completely different behavior arises. We assume that the action of the hunting neuron overrides the effect of the statolith and thus $`S_i0,i`$. The dynamical system (1) with the $`\rho _{ij}`$ values specified above (see also Fig. 7) and with a stimuli from the hunting neuron given, for example, by $`H_i=(0.730,0.123,0.301,0.203,0.458,0.903)`$ has a strange attractor in the phase space (see Fig. 8). This means that the SRN network generates new information (a chaotic signal with positive Kolmogorov-Sinai entropy) in the presence of the prey, which controls the CPGs and, in fact, organizes the random-like behavior of Clione.

The origin of the chaoticity in such dynamical system can be explained in the following manner \[Varona et al., 2001\]: due to the diversity in the strengths of the inhibitory connections we may consider the complete network as two weakly coupled WLC triangle networks. Independently each of them has a closed heteroclinic contour (see Fig. 7), which becomes a limit cycle under the action of a small perturbation (see Sec. 4). The periodic oscillations corresponding to these limit cycles have, in general, different frequencies that are extremely sensitive to the distance to the heteroclinic loop in the non perturbed system (such oscillations are strongly non-synchronous). As we showed the weak interaction of these WLC triangles (nonlinear oscillators) generate chaos in wide regions of the control parameter space. New experiments have confirmed the validity of the model and its predictions \[Levi et al., 2002\].

## 6 Discussion

The stimulus dependent sequential switching of neurons or group of neurons (clusters), named WinnerLess Competition, is able to solve the fundamental contradiction between sensitivity and robustness of the sensory recognition. The key points on which the WLC networks are based are: (i) the heteroclinic contour corresponding to a specific sequence of switching has a large basin of attraction, i.e. a specific sequence is stable; and (ii) the topology of the heteroclinic contour sensitively depends on the incoming signals, i.e. high resolution or sensitivity. Both features are actually present only in the case if under the action of a perturbation the discussed heteroclinic contour, which is structurally unstable, is transformed to the limit cycle in its vicinity with the same topology as the contour. In this paper we have discussed the conditions for such topological stability, and we have showed that computing with separatrices based on the WLC principle is a very natural and powerful strategy for information processing in real neural systems. Any kind of sequential activity can be programmed by a network with stimulus dependent nonsymmetric inhibitory connections. It can be the creation of spatio-temporal patterns of motor activity, the transformation of the spatial information into spatio temporal information for successful recognition and many other computations. In addition, we wish to mention that two important computational functions can be successfully implemented by computation with separatrices. These are: (i) sequential memory storage, and (ii) feature binding.

In reference \[Seliger et al., 2003\] the authors suggest a new biologically-motivated model of sequential spatial memory which is based on the WLC principle. Each stimulus event (visual image, odor, etc…) is represented by a saddle point in the phase space of the system, and a network of one-dimensional separatrices leads the system along the sequence of events in the specific episode. After the learning process, such system is capable of an associative retrieval of the pre-recorded sequence of spatial patterns.

A binding problem occurs when two (or more) different events, e.g. scenes, features, or behaviors are represented by different neural ensembles simultaneously, and for some reason they are all connected with each other. Eventually, these coherent features are integrated by the nervous system of the animal onto a perceptual object, even if the features are dispersed among different sensory systems or subsystems. The binding is ubiquitous and occurs whenever a simultaneous remembrance or representation is important. The most common approach in the modeling of binding is to involve time in operation (von der Malsburg, Singer, and others). The idea is to use the coincidence of certain events in the dynamics of different neural units for binding. This is a dynamic binding. Usually, dynamic binding is represented by synchronous neurons or neurons that are in resonance with and external field. However, dynamical events like phase or frequency variations usually are not very reproducible and robust. It is reasonable to hypothesize that brain circuits that display sequential switching of neural activity \[Abeles et al., 1995\] use the coincidence of this switching to implement dynamic binding of different WLC networks.

In the conclusion we have to emphasize that for large inhibitory neural ensembles it is not necessary to have specific connections that satisfy the conditions formulated above for the existence and stability of the WLC dynamics. If the connections are random, the neurons in the ensemble can be separated in three groups: (i) neurons that are weakly coupled with others (they behave like nearly independent elements), (ii) neurons with strong but close to symmetric connections (they form a Hopfield like network just with simple attractors), and (iii) neurons with nonsymmetric connections that demonstrate the WLC dynamics. Because the WLC dynamics is a structurally stable phenomenon, it is reasonable to hypothesize that the perturbation of the third group by the first two does not destroy the switching activity. Recent computer experiments \[Zhigulin, 2002\] have confirmed this hypothesis.

ACKNOWLEDGMENTS

Support for this work came from NIH grant 2R01 NS38022-05A1, Department of Energy grant DE-FG03-96ER14592 and NSF/EIA-0130708. V.A. was supported by CONACyT grant 485100-3-36445-E and by UC MEXUS-CONACyT grant. P.V. was supported by MCyT BFI2000-0157.


format: text

text:
# The three-dimensional Anderson model of localization with binary random potential

## I Introduction

Amorphous metallic alloys (AMA) offer the possibility of continuously changing their composition while at the same time avoiding structural phase transformations. Thus they allow for systematic studies of their physical properties within a single phase as temperature and other external control parameters, e.g., pressure, are varied. Many such investigations (see, e.g., Has83 ) have been devoted to the investigation of the electrophysical properties (in particular, electrical resistivity) of AMA; these studies have revealed considerable differences in the behavior of the electrical resistivity $`\varrho `$ of AMA in comparison with that of their crystalline counterparts.

The transport properties of transition-metal-based AMA are of special interest because of unique physical anomalies, unexplained by conventional transport theory, that exist in these materials. Among them are the negative temperature coefficient of resistivity, the Mooij correlation of resistivity and thermoelectric power Moo73  — the higher the resistivity the lower is the temperature coefficient of the resistivity —, the resistivity saturation, the sign reversal of the Hall coefficient, a negative magneto-resistance, and the breakdown of Mattiessen’s rule. For a review, we refer the reader to ref. LeeR85 . The first attempt for a theoretical understanding of the electronic transport in noncrystalline materials was the Ziman theory and its extensions Has83 . However, the applicability of the Ziman diffraction model for strongly disordered and thus high-resistivity amorphous alloys is questionable CotM82 ; MeiC84 . In addition, the Ziman theory cannot explain, e.g., that the changes of the resistivity $`(\varrho _{\mathrm{liquid}}\varrho _{\mathrm{solid}})/\varrho _{\mathrm{solid}}`$ at the melting point are only $`0.01`$ and $`0.09`$ for Fe and Co, respectively Has83 . Thus it appears that the resistivity of these two elements depends only weakly on their atomic structure factor $`S(k)`$.

At low temperature $`T`$, an even more significant difference between the behavior of crystals on the one hand and disordered solids on the other hand is seen: sufficiently strong disorder can give rise to a transition of the transport properties from conducting behavior with resistance $`R>0`$ to insulating behavior with $`R=\mathrm{}`$ as $`T0`$ as was pointed out by Anderson in 1958 And58 , This phenomena is called the disorder-driven metal-insulator transition (MIT) BelK94 ; KraM93 ; LeeR85  and it is characteristic to non-crystalline solids. The mechanism underlying this MIT was attributed by Anderson not to a finite gap in the energy spectrum which is responsible for an MIT in band gap or Mott insulators AshM76 . Rather, he argued that the disorder will lead to interference of the electronic wave function $`\psi (𝐫)`$ with itself such that it is no longer extended over the whole solid but is instead confined to a small part of the solid. This localization effect excludes the possibility of diffusion at $`T=0`$ so that the system is an insulator. A highly successful theoretical approach to this disorder-induced MIT was put forward in 1979 by Abrahams et al. AbrALR79 . This “scaling hypothesis of localization” details the existence of an MIT for non-interacting electrons in three-dimensional disordered systems at zero magnetic field and in the absence of spin-orbit coupling.

These ideas have also been applied to the analysis of AMA transport properties Imr79 . They provide a background for an explanation of the Mooij correlation ParSK01  and were successful in the description of the transport properties of amorphous semiconductors Mad78 . Many papers have since discussed the influence of localization (or quantum interference effects) on the transport properties of AMA HicGH86 ; BogH90 ; HowHM88 ; BieFGO84 . It was shown SeeLHM99  that Anderson localization is responsible for regions of high resistivity. The importance of quantum corrections in analyzing the conductivity of so-called highly resistive ($`\varrho >150\mu \mathrm{\Omega }`$cm) metallic glasses was demonstrated HowG86 . There are two sources of these corrections: the disorder-induced “localization effect” and the electronic “interaction effect”. The former caused the observed higher temperature resistivity. Furthermore, it has been argued RosHPK97 ; Tsu86  that for amorphous metals localization effects are valid even at room temperature.

In a recent study NakPSZ00  the influence of disorder on the transport properties of Ti-Ni and Ti-Cu metallic glasses was analyzed within the Anderson model of localization. The results elucidated qualitatively the correlation between the observed details of the electronic structure and the behavior of resistivity vs. temperature for this binary AMA. It was shown that the temperature coefficient of resistivity depends on the position of the Fermi level $`E_\mathrm{F}`$ and on the (assumed) position of the mobility edge $`E_\mathrm{c}`$. The latter separates extended (conducting) states from localized (insulating) states KraM93 ; GruS95 . Unfortunately, the value of $`E_\mathrm{c}`$ had to be inferred from the hypothesis that the two subbands of the valence band of the binary alloy have tails with localized electronic states like the tails of the usual one-band Anderson model And58 . Evidently, a calculation of the position of the mobility edge $`E_\mathrm{c}`$ should be carried out for the case of two subbands. This is what we intend to do in the present manuscript. In particular, we will investigate whether localized electronic states exist in the central region between the two subbands as assumed in NakPSZ00 . Furthermore, we show how the concentration parameter, absent in the usual single-band Anderson model, influences the position of the mobility edges.

## II Numerical approach

### II.1 The binary Anderson model of localization

According to the photo-emission data and theoretical estimates Haf91 ; NakPSZ00 , the valence band of a binary, transition-metal-based amorphous alloy can be assumed to be a superposition of two valence bands $`A`$ and $`B`$. In order to model this within the Anderson model, we use the standard Anderson Hamiltonian And58

$$𝐇=\underset{i}{}\epsilon _i|ii|+\underset{ij}{}t_{ij}|ij|$$
(1)
with orthonormal states $`|i`$ corresponding to electrons located at sites $`i=(x,y,z)`$ of a regular cubic lattice with periodic boundary conditions. The hopping integrals $`t_{ij}`$ are non-zero only for nearest neighbors and we set the energy scale by choosing $`t_{ij}=t=1`$. The two-subband structure observed in the experiments is incorporated into the potential energies $`\epsilon _i`$. The randomness is modeled by

1. independent random variations $`\epsilon _{i(A)}[\epsilon _AW_A/2,\epsilon _A+W_A/2]`$ and $`\epsilon _{i(B)}[\epsilon _BW_B/2,\epsilon _B+W_B/2]`$,
2. random spatial distribution of potential energies $`\epsilon _A`$ and $`\epsilon _B`$ within the cubic lattice.

The mean values $`\epsilon _A`$ and $`\epsilon _B`$ are chosen according to the central energies of the two subbands. The parameters $`W_A`$ and $`W_B`$ specify the disorder strength in each energy band and the random variation of $`A`$ and $`B`$ sites models the compositional disorder of the binary alloy. Note that with this choice of parameters, the band edges of the subband $`A`$ in the limit of large system size are given by $`\epsilon _A6W_A/2`$, $`\epsilon _A+6+W_A/2`$ and similarly for the other subband with $`AB`$.

In summary, the model is described by the following parameters

1. $`\epsilon _{AB}=\epsilon _B\epsilon _A`$, the distance between the centers of the two bands (we choose $`\epsilon _B>\epsilon _A`$ so that the $`A`$ band is the lower one),
2. $`x_A`$, the concentration of $`A`$ sites (of course, then $`x_B=1x_A`$ is the concentration of the $`B`$ sites),
3. $`W_A`$ and $`W_B`$, the widths of the box distributions of the potential energy.

### II.2 Transfer-matrix method

Since we are interested in the position of the mobility edges and thus the localization lengths, we use the transfer-matrix method (TMM) KraM93 ; PicS81a ; MacK81  for the investigation of model (1). The localization length $`\lambda `$ describes the exponential decay of the wave function and we compute it using TMM for quasi-1D bars of cross section $`M\times M`$ and length $`LM`$. As is customary in the TMM approach, the stationary Schrödinger equation $`𝐇\mathrm{\Psi }=E\mathrm{\Psi }`$ is rewritten in the recursive form

$$\left(\genfrac{}{}{0pt}{}{\mathrm{\Psi }_{i+1}}{\mathrm{\Psi }_i}\right)=\left(\begin{array}{cc}E\mathrm{𝟏}𝐇_i& \mathrm{𝟏}\\ \mathrm{𝟏}& \mathrm{𝟎}\end{array}\right)\left(\genfrac{}{}{0pt}{}{\mathrm{\Psi }_i}{\mathrm{\Psi }_{i1}}\right)=𝐓_i\left(\genfrac{}{}{0pt}{}{\mathrm{\Psi }_i}{\mathrm{\Psi }_{i1}}\right).$$
(2)
$`\mathrm{\Psi }_i`$, $`𝐇_i`$, and $`𝐓_i`$ are wave function, Hamiltonian matrix, and transfer matrix of the $`i`$th slice of the bar, respectively. Unit and zero matrices are denoted by $`\mathrm{𝟏}`$ and $`\mathrm{𝟎}`$. Given an initial condition $`\left(\genfrac{}{}{0pt}{}{\mathrm{\Psi }_1}{\mathrm{\Psi }_0}\right)`$ equation (2) allows a recursive construction of the wave function in the bar geometry by adding more and more slices. $`\lambda (M,w)`$ is then obtained from the smallest Lyapunov exponent of the product $`𝐓_L𝐓_{L1}\mathrm{}𝐓_2𝐓_1`$ of transfer matrices MacK83 , where the length $`L`$ of the bar is increased until the desired accuracy of $`\lambda `$ is achieved. With increasing cross section of the bar the reduced localization length $`\mathrm{\Lambda }_M(w)=\lambda (M,w)/M`$ decreases for localized states and increases for extended states. Thus it is possible to determine the critical parameter $`w_\mathrm{c}`$ at which $`\mathrm{\Lambda }_M`$ is constant as a function of the varied parameter $`w`$ — e.g., $`E`$, $`\epsilon _{AB}`$, $`x_A`$, $`W_A`$, or $`W_B`$ — from plots of $`\mathrm{\Lambda }_M`$ versus $`M`$.

### II.3 Finite-size scaling

The MIT in the Anderson model of localization is expected to be a second-order phase transition BelK94 ; AbrALR79 . It is characterized by a divergent correlation length $`\xi _{\mathrm{}}(w)=C|ww_\mathrm{c}|^\nu `$, where $`\nu `$ is the critical exponent, $`C`$ is a constant KraM93 , $`w`$ is any of the external control parameters given above and $`w_\mathrm{c}`$ is its critical value at which the MIT occurs. To construct the correlation length of the infinite system $`\xi _{\mathrm{}}`$ from finite-size data $`\mathrm{\Lambda }_M`$ ZamLES96a ; KraM93 ; PicS81a ; MacK81 , the one-parameter scaling hypothesis Tho74  is employed,

$$\mathrm{\Lambda }_M=f(M/\xi _{\mathrm{}}).$$
(3)
All values of $`\mathrm{\Lambda }_M(w)`$ are expected to collapse onto a single scaling curve $`f`$, when the system size is scaled by $`\xi _{\mathrm{}}(w)`$. In a system with MIT such a scaling curve consists of two branches corresponding to the localized and the extended phase. One might determine $`\nu `$ from fitting $`\xi _{\mathrm{}}`$ by a finite-size scaling (FSS) procedure MacK83 . But a higher accuracy can be achieved by fitting directly the raw data MacK83 . We use fit functions SleO99a  which include two kinds of corrections to scaling: (i) nonlinearities of the $`w`$ dependence of the scaling variable and (ii) an irrelevant scaling variable with exponent $`y`$. Specifically, we fit

$$\mathrm{\Lambda }_M=\stackrel{~}{f}_0(\chi _\mathrm{r}M^{1/\nu })+M^y\stackrel{~}{f}_1(\chi _\mathrm{r}M^{1/\nu }),$$
(4)
$$\stackrel{~}{f}_n=\underset{i=0}{\overset{n_\mathrm{r}}{}}a_{ni}\chi _\mathrm{r}^iM^{i/\nu },\chi _\mathrm{r}(\omega )=\omega +\underset{n=2}{\overset{m_\mathrm{r}}{}}b_n\omega ^n$$
(5)
with $`\omega =|w_\mathrm{c}w|/w_\mathrm{c}`$ and expansion coefficients $`a_{ni}`$ and $`b_n`$. Choosing the orders $`n_\mathrm{r}`$ and $`m_\mathrm{r}`$ of the expansions larger than one, terms with higher order than linear in the $`w`$ dependence appear. This allows to fit a wider $`w`$ range around $`w_\mathrm{c}`$ than with the previously used linear fitting Mac94 . The linear region is usually very small. The second term in equation (4) describes the systematic shift of the crossing point of the $`\mathrm{\Lambda }_M(w)`$ curves Mac94 ; SleO99a .

In the present case of the two-band Anderson model, we can in principle have an MIT as a function of

1. energy $`E`$ for fixed $`W_A`$, $`W_B`$, $`\epsilon _{AB}`$ and $`x_A`$,
2. disorder strengths $`W_A`$, $`W_B`$ for fixed $`E`$, $`\epsilon _{AB}`$ and $`x_A`$,
3. concentration $`x_A`$ for fixed $`E`$, $`\epsilon _{AB}`$ and $`W_A`$, $`W_B`$.

Due to universality, we expect the corresponding critical exponents $`\nu _E`$, $`\nu _{W_A}`$, $`\nu _{W_B}`$ and $`\nu _x`$ to be the same. Additionally for each control parameter we test whether the fitted values of $`w_\mathrm{c}`$ and $`\nu `$ are compatible when using different expansions of the fit function, i.e., different orders $`n_\mathrm{r}`$ and $`m_\mathrm{r}`$ MilRS00 .

## III Results and discussion

For the determination of the mobility edge GruS95 , we have determined the localization lengths from TMM with an error equal to or less than $`1\%`$ for cross-sections up to $`M=12`$. For the computation of the critical exponents, we have used the available $`1\%`$ data and generated additional data with $`0.05\%`$ error up to at most $`0.1\%`$ at selected points in the phase diagrams for high-precision estimates.

In Fig. 1, we show a typical dependence of the reduced localization length $`\mathrm{\Lambda }_M`$ on energy for different system sizes $`M`$. The two peaks correspond roughly to the positions of the two subbands at $`\epsilon _A`$ and $`\epsilon _B`$, respectively. This is in accordance with the theoretical and experimental results NakPSZ00 ; Haf91  mentioned in the introduction as motivation of our studies. Note that for the chosen parameters, the density-of-states limit given in section II.1 for the highest energy state in the lower subband coincides with the smallest possible energy state in the upper subband at $`\epsilon _A+6+W_A/2=\epsilon _B6W_B/2`$.

At the positions indicated by the dashed lines in Fig. 1, we see a reversal in the systematic size dependence of the reduced localization lengths. Therefore these indicate the mobility edges and we observe localized electronic states at the band edges and in the central region between the $`A`$ and $`B`$ subbands. We remark that the $`\mathrm{\Lambda }_M(E)`$ curves demonstrate a peculiar feature: the positions of the maximum values of $`\mathrm{\Lambda }_M`$ do not coincide with the central energies $`\epsilon _A`$ and $`\epsilon _B`$ of the two subbands. Thus although the density of states between the two bands is appreciable, we observe the surprising fact that the region of extended states in each band has shrunk.

In Fig. 2 we show high precision data at the upper mobility edge of the $`A`$ subband at $`E_\mathrm{c}2`$. Again the accuracy of the data is so high that error bars would be smaller than the symbols.

Data and FSS analysis of similar quality will be used when estimating critical exponents in the following.

### III.1 Energy-disorder phase diagram

Let us now investigate how the positions of the mobility edges change when $`W_A`$ and $`W_B`$ are changed. We set $`\epsilon _{AB}=17`$ and $`x_A=x_B=0.5`$. For convenience we choose $`W_AW_B`$. This leads to a symmetry for the energy dependences of $`\mathrm{\Lambda }_M(E)`$ between lower and upper subband as shown in Fig. 1. Consequently, the phase diagram is symmetric with respect to $`\epsilon _A+\epsilon _{AB}/2`$. We find extended electronic states in the vicinity of $`\epsilon _A`$ and for $`W_AW_\mathrm{c}6.5`$ (instead of $`16.5`$ as for the usual Anderson model) as presented in Fig. 3. For larger $`W_A`$ or larger $`|E\epsilon _A|`$, the states are localized.

Next, we keep $`W_A=4`$ fixed and vary $`W_B`$ from 0 to 8. We find for $`\epsilon _{AB}=17`$ that the position of mobility edges in the subband $`A`$ is only slightly modified. The same behavior is observed for $`\epsilon _{AB}=8`$, when the overlap of the two subbands is large.

### III.2 Energy-$`\epsilon _{AB}`$ phase diagram

Keeping $`x_A=0.5`$ and fixed disorder $`W_A=W_B=4`$, we now vary the subband spacing $`\epsilon _{AB}`$. From an experimental point of view, this is the difference between the atomic constituents of the alloy or, more precisely, between their ionization energies. For a variety of $`\epsilon _{AB}`$ values we have obtained $`\mathrm{\Lambda }_M`$ curves which are qualitatively very similar to Fig. 1. However, for $`\epsilon _{AB}6`$, no localized electronic states in the central region can be found. In Fig. 4, we show the mobility edges for different $`\epsilon _{AB}`$.

We note that the mobility edge for the usual Anderson model equals approximately $`\pm 6.2`$ at disorder $`W=4`$ KraM93 ; BulSK87 . This agrees with the result for $`\epsilon _{AB}=0`$ in Fig. 4.

### III.3 Energy-concentration phase diagram

Obviously, the values for $`\mathrm{\Lambda }_M`$ and the mobility edges should strongly depend on the concentration. This is indeed the case as we show in Fig. 5. Already a $`10\%`$ difference in atomic composition leads to a pronounced asymmetry of the $`\mathrm{\Lambda }_M`$ curves (compare, e.g., with Fig. 1).

The energy-concentration phase diagram for $`\epsilon _{AB}=8`$, and $`W_A=W_B=4`$ is presented in Fig. 6. For $`x_A<0.32`$ ($`x_B>0.68`$), all electronic states of subband $`A`$ ($`B`$) are localized. Figs. 5 and 6 correspond to the case of overlapping bands such that the center of subband $`A`$ coincides with the edge of subband $`B`$ and vice versa.

### III.4 Concentration-$`\epsilon _{AB}`$ phase diagram

The concentration-$`\epsilon _{AB}`$ phase diagram is presented in Fig. 7 for $`E=\epsilon _A=0`$ and $`W_A=W_B=4`$.

We see that the critical concentration is larger than the classical percolation threshold $`x_A=0.3116`$ for large $`\epsilon _{AB}`$. We interpret this as being due to the localization effect in this quantum situation: a larger cluster of connected sites is needed to support extended states on sites of the sublattice $`A`$ (or $`B`$). For very large $`\epsilon _{AB}=100`$ and $`10^3`$ (not shown) the critical concentration decreases to $`0.416`$ and $`0.384`$, respectively, but remains larger than in the classical case. For large subband spacing sites with onsite potential $`\epsilon _B`$ can be viewed as potential barriers to electrons from subband $`A`$ and vice versa. Thus electrons with energy $`\epsilon _A`$ have a low probability to hop onto $`B`$ sites and we are effectively studying the case of percolation. But even though $`\epsilon _{AB}`$ might be large, there will always remain a finite but small probability for tunneling thus distinguishing our model from classical percolation. And since we are studying the wave equation (1), localization effects not present in classical percolation are important and lead to insulating behavior even when the concentration of sites supports a classically percolating cluster. Therefore, we find that the threshold for transport (extended states) in the quantum case is higher than for the classical one. We note that similar effects have been observed previously in studies of quantum percolation AviL92 .

On the other hand, for small $`\epsilon _{AB}`$, the critical concentration becomes less than the percolation threshold since tunneling between $`A`$ and $`B`$ sites is more significant. The phase diagram in Fig. 7 does not exactly show the mobility edges at which all extended states disappear for a given $`\epsilon _{AB}`$, because it is determined at $`E=0`$. But as shown in Figs. 1 and 6, the $`A`$ band is not symmetric with respect to $`\epsilon _A(=0)`$. As a consequence, more extended states occur for energies just below $`\epsilon _A`$. This means that the phase boundary in Fig. 7 will be shifted slightly downwards, if $`E`$ is varied, too.

We note that an analogous phase diagram has been obtained recently NazBR02  for the case of a two-dimensional binary alloy. The shape of the curve separating localized and extended electronic states is quite similar to our result shown in Fig. 7.

### III.5 Disorder-concentration phase diagram

The concentration dependence of the critical disorder — the disorder at which all states for a given energy become localized — is presented in Fig. 8 for energy $`E=\epsilon _A`$. One can see that the critical disorder strongly depends on the concentration. For $`x_A=1`$, we recover the result $`W_\mathrm{c}=16.5`$ KraM93  of the single-band Anderson model.

No extended electronic states are observed for $`x_A<0.35`$ in agreement with the idea that for small concentration the electronic states of these atoms form the usual localized donor subband of separated impurities.

For the concentration interval from $`0.35`$ to $`0.40`$ we obtain an additional transition as a function of disorder at small $`W_A`$ besides the usual Anderson MIT. For example, as shown in Fig. 9 for $`x_A=0.38`$ we observe localized behavior of $`\mathrm{\Lambda }_M(M)`$ for small disorders. Increasing the disorder beyond $`W_A>0.7`$, we see extended states. At $`W_A3.2`$, the character of the states changes back to localized.

This behavior can be understood as follows: for $`x_A0.38`$ as indicated by the arrow in Fig. 8, the system at $`W_A=0`$ represents a binary alloy with large energy separation $`\epsilon _{AB}=17`$. All states are localized, i.e., non-quantum percolating, on an already classically percolating set of sites $`\{i(A)\}`$ in the $`A`$ band. Small additional potential disorder will lead to an increase of certain $`\epsilon _{i(A)}`$ as well as a decrease in certain $`\epsilon _{i(B)}`$. A small change of the scattering conditions in these sites can then quickly allow enhanced transport across a now quantum-percolating backbone. But upon further increasing the disorder, the localization in each band will quickly lead to localization again.

### III.6 Value of the critical exponent

The critical exponent for the present model can be estimated by crossing the mobility edges as function of either disorder, energy or concentration, giving rise to the three exponents $`\nu _W`$, $`\nu _E`$ and $`\nu _x`$, respectively. Due to universality, their values should coincide with the critical exponent $`\nu 1.6`$ of the single-band Anderson model since the universality class remains unchanged by the introduction of the binary disorder HofS93b .

We have estimated the critical exponent $`\nu `$ for almost all points in the phase diagrams 3, 4, 6, 7 and 8 as described in section II.3. For example, we compute $`\nu _W`$ along the mobility edge of the disorder-concentration phase diagram of Fig. 8. In Fig. 10 we show the resulting estimates when using different orders $`n_r`$ and $`m_r`$ for the fit function (4). As usual, the spread in values is somewhat larger than the least-square-error bars seem to suggest MilRS00 . Summarizing the results, we find $`\nu _W=1.56\pm 0.06`$, $`\nu _E=1.64\pm 0.05`$ and $`\nu _x=1.60\pm 0.07`$, compatible with $`\nu 1.6`$ MilRSU00 .

## IV Discussion and comparison with experiments

Let us now compare the above results of the binary Anderson model with experimental measurements of transport properties of AMA NakPSZ00 ; Cot76 . Such a comparison will of course be purely qualitative. Nevertheless, it already suffices to understand much of the transport properties of such AMA. E.g., as mentioned in section I, it is the position of the mobility edge relative to the value of the Fermi energy that determines the properties of transition-metal-based AMA. In this spirit, the Mooji correlation Moo73  is due to $`E_\mathrm{F}`$ being situated in a region of localized states, whereas the weak changes in resistivity at melting for Fe and Co as well as the positive temperature coefficient of resistivity for Fe-based AMA result from $`E_\mathrm{F}`$ being located in the center of the 3d-electronic band and hence, in the region of extended states.

### IV.1 Transport properties of Ni<sub>1-x</sub>P<sub>x</sub>

The ionization energy of P is higher than that of Ni, so our $`A`$ band corresponds to P atoms and the $`B`$ band to Ni atoms with $`\epsilon _B>\epsilon _A`$ as, e.g., shown in Fig. 6. According to the theoretical and experimental results of Ref. BabPSZ95 , $`E_\mathrm{F}`$ of the amorphous alloy Ni<sub>1-x</sub>P<sub>x</sub> is located near the upper edge of the Ni-3d band (the $`B`$ band in our study). When the P concentration $`x_\mathrm{P}`$ increases from $`0.15`$ to $`0.27`$ the temperature coefficient of resistivity (TCR) decreases from positive to negative values Cot76 . The MIT is observed near $`x_\mathrm{P}=0.25`$. In Fig. 6, we indicated the corresponding region in the phase diagram. The two-band, binary random potential Anderson model can indeed qualitatively capture the observed experimental results. Even a quantitative comparison might be possible if a reliable way of obtaining values for $`\epsilon _A`$, $`\epsilon _B`$, and $`W_A`$ and $`W_B`$ becomes known. Note, however, that in the experiments, it is the hopping between $`A`$ and $`B`$ atoms at their respective random positions that leads to the broadening of the $`A`$ and $`B`$ bands and we have modeled this only qualitatively in the present work by the onsite disorders $`W_A`$ and $`W_B`$.

### IV.2 Transport properties of Ti<sub>1-x-y</sub>Cu<sub>x</sub>A<sub>y</sub> and Ti<sub>1-x-y</sub>Ni<sub>x</sub>A<sub>y</sub>

A previous experimental study of the electronic structure of Ti<sub>1-x-y</sub>Cu<sub>x</sub>A<sub>y</sub> and Ti<sub>1-x-y</sub>Ni<sub>x</sub>A<sub>y</sub> amorphous alloys NakPSZ00  (the symbol A<sub>y</sub> represents additional admixtures) has shown that the valence band has two main peaks. The lower peak is formed predominantly by the 3d states of Cu or Ni whereas the upper peak is due to the 3d states of Ti. The Fermi level is located in the central region between the peaks and is hardly shifted upon changing the concentrations.

Within the present two-band model, this experimental situation may be modeled as shown in Fig. 6. The line for the Ti<sub>1-x-y</sub>Ni<sub>x</sub>A<sub>y</sub> alloys is situated somewhat below the one for Ti<sub>1-x-y</sub>Cu<sub>x</sub>A<sub>y</sub> because Ni has one electron less than Cu. We believe that it is this difference that leads to the smaller value of the negative TCR for Ti<sub>1-x-y</sub>Ni<sub>x</sub>A<sub>y</sub> when compared to Ti<sub>1-x-y</sub>Cu<sub>x</sub>A<sub>y</sub>. In both cases, $`E_\mathrm{F}`$ lies in a region of localized states in agreement with the negative TCR. In Table 1, we show the values obtained for TCR.

Let us now analyze the concentration dependence of the TCR. When the Ti concentration decreases, the TCR also decreases as shown in Table 1. This behavior can be explained from the concentration-energy phase diagram of Fig. 6 by the increase of the mobility edge $`E_\mathrm{c}`$. This in turn leads to an increasing distance between $`E_\mathrm{F}`$ and the region of unoccupied extended electronic states. The states close to $`E_\mathrm{F}`$ become more localized, leading to a more negative TCR.

Therefore, at least at a qualitative level, the present two-band Anderson model explains not only the negative TCR of Ti<sub>1-x-y</sub>Cu<sub>x</sub>A<sub>y</sub> and Ti<sub>1-x-y</sub>Ni<sub>x</sub>A<sub>y</sub> glasses but also the changes in TCR due to changes in composition (Cu or Ni) and concentration. Furthermore, metallic glasses with a “metal”-like conductivity can be also treated within this model. For example, the Fermi level of Fe<sub>86</sub>B<sub>14</sub> alloy is situated in the region of extended states as indicated in Fig. 6 in agreement with its positive TCR NakPSZ00 .

### IV.3 Overlapping the two bands

Thus far, motivated by the experimental results, we concentrated on the case of separated bands, i.e., large $`\epsilon _{AB}`$. To make the two regions of extended states in Fig. 6 overlap, we can in principle (i) keep $`\epsilon _{AB}`$ fixed and change the disorder, or (ii) keep fixed disorder and decrease $`\epsilon _{AB}`$ as in Fig. 4. From Fig. 3, we know that the position of mobility edges does not change a lot when decreasing $`W_A`$. On the other hand, increasing $`W_A`$ will lead to a narrower region of extended states in Fig. 6 and no overlap. Therefore, we choose to reduce $`\epsilon _{AB}`$ in order to study what happens when the two bands overlap in the energy-concentration diagram.

In Fig. 11, we show such an overlap for $`\epsilon _{AB}=6`$ and $`W_A=W_B=4`$.

This situation corresponds to large overlap of the two subbands. We see that the shape of the regions of extended states remains nearly unchanged. This leads to non-monotonicity in the behavior of the upper and lower mobility edges. We remark that in order to resolve this behavior, the accuracy of the TMM had to be increased significantly up to very small errors of $`0.05\%`$ for the localization lengths. We emphasize that a similar accuracy was needed to resolve the non-monotonicity in the disorder-concentration diagram in Fig. 8.

When further decreasing $`\epsilon _{AB}`$ to $`5`$, the behavior of the mobility edges becomes regular again and their values eventually reduce to the values of the well-known phase diagram of the single-band Anderson model BulKM85 .

## V Conclusions

We have shown that the three-dimensional Anderson model of localization with binary random potential disorder allows to explain not only various peculiarities of the transition-metal-based amorphous alloys but also how these change as the chemical composition changes. Even for a comparatively large overlap of the two subbands, a region of localized electronic states exists in the central part between the two subbands. This confirms the assumption made in Ref. NakPSZ00  for the experimental results in AMA. Even a quantitative analysis, using suitably extracted parameters from density-of-states measurements appears possible. Of course in this case ternary and quaternary AMA will have to be distinguished by additional onsite potentials similar to $`\epsilon _A`$ and $`\epsilon _B`$ for the binary model.

In addition to the usefulness of the two-band Anderson model for a qualitative comparison with experiments, there is also the interesting fact that the model allows to study the concentration dependence and the associated MIT. This concentration dependence is clearly a prominent feature present in most disordered materials. We find that, as expected from universality, the universality class of this additional transition remains unchanged with respect to a single-band Anderson model.

###### Acknowledgements.

We thank A. Eilmes, S. Klassert and P. Stollmann for useful comments. This work was supported by grants from the DAAD via Ref. 322, A/01/20012 and the DFG via SFB393 .