text:
# Doubled Full Shot Noise in Quantum Coherent Superconductor - Semiconductor Junctions

Present adress \] Max Planck Institut für Festkörper-
forschung, Heisenbergstr. 1, D-70569 Stuttgart, Germany.

## Abstract

We performed low temperature shot noise measurements in superconductor (TiN) - strongly disordered normal metal (heavily doped Si) weakly transparent junctions. We show that the conductance has a maximum due to coherent multiple Andreev reflections at low energy and that the shot noise is then twice the Poisson noise ($`S=4eI`$). When the subgap conductance reaches its minimum at finite voltage the shot noise changes to the normal value ($`S=2eI`$) due to a large quasiparticle contribution.

We know, from early measurements, that the full shot noise in an electronic device $`S_{Poisson}=2qI`$ (first measured by W. Schottky in a vacuum diode) is proportional to the mean value of the current $`I`$ and to the charge $`q`$ of the carriers Schottky . This result holds for N-I-N tunnel junctions where N is a normal metal and I an insulating barrier Blanter00  with $`q=e`$ the electronic charge. In S-I-N, due to electron pairing in the superconductor (S), the shot noise is expected to be twice the full shot noise : $`S=4eI`$. However, in such junctions the subgap current is very small and the shot noise is not measurable. The subgap current can be restored if the quasiparticles of the normal metal are coherently backscattered towards the interface. This reflectionless tunneling regime can be achieved by adding a second barrier in the normal part of the junction (S-I-N-I-N) or when the normal metal is disordered enough Quirion01 ; Kastalsky91 ; VanWees92 ; Quirion00 ; Beenakker97 ; Volkov94 ; Bakker ; Magnee ; Beltram . The enhancement of the subgap current is only seen at low energy when the electron-hole coherence time in the normal metal is longer than the time it takes for a quasiparticles to return to the interface. Then, the coherent addition of two (or more) Andreev reflections, each of them with a very small probability $`\mathrm{\Gamma }^2`$ ($`\mathrm{\Gamma }`$ is the transparency of the barrier), yields to an increase of the Andreev current through the interface. This effect can be large, leading to an Andreev current proportional to $`\mathrm{\Gamma }`$ instead of $`\mathrm{\Gamma }^2\mathrm{\Gamma }1`$ and can be comparable to the normal current above the gap (also proportional to $`\mathrm{\Gamma }`$). Another way to increase the subgap current is to use highly transparent S-N junctions. In this case doubled shot noise is predicted and has been observed experimentally at various temperatures Kozhevnikov00 ; Jehl00 . However, the noise level was always three times smaller than the doubled full shot noise because of the diffusive nature of the normal metal used in these experiments. Moreover, as shown experimentally and reproduced theoretically Nagaev01 , the doubling of the shot noise occurs at any energies below the superconducting gap and electron-hole coherence is not required.

In this letter, we report shot noise measurements in a junction where a superconductor (TiN) is in contact with heavily doped silicon. We show that the shot noise is twice the full shot noise at low energy and equals the Poisson value at bias much smaller than the superconducting gap. This behavior evidences a crossover from a low bias Andreev dominated to a large bias quasiparticle dominated subgap conductance.

The sample is made of two 100/10 nm thick TiN/TiSi2 contacts on top of a silicon substrate. The silicon is heavily doped ($`n_e=2.10^{19}cm^3`$) over a thickness of $`0.6\mu m`$ Quirion00 . The two contacts are squares of 1 mm wide and the distance between them is $`1\mu m`$ or $`2\mu m`$. We present the results obtained with the $`1\mu m`$ sample and found the same results for the $`2\mu m`$ sample. Because of the doping, the silicon is metallic ($`k_Fl_e3`$). Moreover, the superconducting TiN/TiSi2 bilayer and the silicon are separated by a Schottky barrier which stays symmetric (non rectifying) in our voltage range. Therefore, the contact is described by a S-I-N junction where I stands for a tunnel barrier. The transport properties have been studied in details by Quirion at al. Quirion00 , and it was shown that the junction presents reflectionless tunneling behavior at low energy ($`T250mK`$ and $`V20\mu V`$). Using reference Volkov94  for fitting the temperature dependence of the zero bias resistance, the following parameters have been obtained: the superconducting gap is $`\mathrm{\Delta }`$ =0.22meV, the damping factor $`\mathrm{\Gamma }_S`$ Dynes78  and the depairing rate $`\gamma _{in}`$ are both relatively large: $`\mathrm{\Gamma }_S/\mathrm{\Delta }=0.15\pm 0.01`$ and $`\gamma _{in}/\mathrm{\Delta }=0.27\pm 0.05`$. The typical Schottky barrier transparency is $`\mathrm{\Gamma }=\mathrm{3.5\; 10}^2`$. The fitting of the subgap conductance implies many parameters and the quantitative agreement should be taken cautiously. As in previous reports of reflectionless tunneling Kastalsky91 ; Bakker ; Magnee ; Beltram  the large background subgap conductance cannot be attributed to Andreev reflection only, but to a large quasiparticule contribution due to a large $`\mathrm{\Gamma }_S`$. Moreover most authors have concluded from conductance measurements to a non uniform barrier interface, that makes the theoretical comparison mostly qualitative. In this context noise measurements are necessary to discriminate between Andreev and quasiparticule contributions.

We also know Quirion00  that coherent effects take place underneath the contacts and do not extend sidewise in the silicon between the contacts. The total sample can therefore be treated as two independent S-I-N junctions connected by a small piece of doped silicon of resistance $`R_{Si}`$. The sample resistance is :

$$R_{sample}=2R_{contact}+R_{Si}=2R_{contact}+N_{sq}R_{sq,Si}$$
(1)
where $`R_{sq,Si}=24\mathrm{\Omega }`$ is the sheet resistance of the silicon and $`N_{sq}=1/1000`$ the number of squares between the two contacts. We also know that the current partially flows below the contacts before entering the superconductor. Therefore, $`R_{contact}`$ ($`0.4\mathrm{\Omega }`$ at low temperature) includes the resistance of the barrier and the resistance of the doped silicon underneath the superconductor which can be larger than in the native film.

In figure 1, we plotted the differential resistance of one contact as a function of the DC current and of the DC voltage drop at the contact, measured with the experimental set up used for the noise. We recover previous results which show that coherent effect appears at energies below $`20\mu eV`$, where the differential resistance shows a maximum Quirion00 . As long as the voltage stays above $`20\mu V`$, the differential resistance increases with decreasing energy because both the quasiparticles and the uncoherent Andreev contributions to the subgap current decrease BTK82 . For low DC voltage and low temperature, the coherence is established between successive reflections and the Andreev current increases. From the temperature and voltage dependences of the conductance it is possible to evaluate the respective parts of the quasiparticle, uncoherent and coherent Andreev contributions. At voltage much larger than $`20\mu V`$, both the Green’s functionsVolkov94  and the BTKBTK82  descriptions give a dominant quasiparticle contribution (the uncoherent Andreev contribution within the BTK model is less than 10 percents). This is due to the large damping factor $`\mathrm{\Gamma }_S`$ for the quasiparticles density of states and to an elevated effective electron temperature $`T_{eff}`$ induced at finite voltageQuirion00 . Such phenomenological $`T_{eff}`$ is introduced to explain why $`eV_{max}k_BT_{max}`$ where $`V_{max}20\mu V`$ is the voltage at which the resistance is maximum and $`T_{max}250mK`$ is the temperature at which the zero bias resistance is maximum. From the theory, such relation is not obeyed if one supposes that the electron temperature is the base phonon temperature. At voltage smaller than $`20\mu V`$, the Andreev contribution (treated within the Green’s functions formalism Volkov94  which describes the coherent part) becomes larger than the quasiparticle part (below $`V=20\mu V`$ the quasiparticle contribution is less than 10 percents). At $`V40\mu V`$, due to an increase of $`T_{eff}`$, both contributions are estimated of the same order.

For shot noise measurements, we used a SQUID based experimental set up Jehl99  which measures the noise of a known macroscopic resistor ($`R_{ref}=0.13\mathrm{\Omega }`$) in series with the sample. The total current noise power $`S_{total}`$ measured, is then given by :

$$S_{total}=2S_{contact}\frac{R_{contact}^2}{\left(R\right)^2}+S_{Si}\frac{R_{Si}^2}{\left(R\right)^2}+S_{ref}\frac{R_{ref}^2}{\left(R\right)^2}$$
(2)
Here, $`S_{contact}`$ is the noise we want to study, $`S_{Si}`$ the thermal noise generated by the silicon between the two contacts, $`S_{ref}`$ the thermal noise of the reference resistor and $`R`$ is the sum of the resistances. We have chosen the sample to be such that $`R_{contact}R_{Si},R_{ref}`$, thus the total noise is principally the noise at the contact and the measurement is mainly insensitive to both the noise between the contacts and to the reference resistor noise. It is also insensitive to a change of the electronic temperature within these parts. At equilibrium, the noise is due to thermal fluctuations ($`S=4k_BT/R`$) and the total measured noise is $`S=4k_BT/(R)`$. Far from equilibrium, equation (2) applies if $`R_{contact}`$ is the differential resistance.

In figure 2, we displayed the current noise power as a function of the DC current that passes through the contact. If we focus on the results obtained at $`T=50mK`$, we clearly see that the shot noise exhibits a kink around $`30\mu A`$ which corresponds to the value of current where the differential resistance shows a broad maximum (figure 1). In the low bias regime ($`I<30\mu A`$), the shot noise is proportional to the current with a slope $`4e`$ corresponding exactly to twice the full shot noise: this is our main result. Note that we do not see the thermal cross-over at $`eVk_BT5\mu eV`$ and that the solid line is $`S=4eI+4k_BT/R_{contact}(V=0)`$. At higher currents ($`I>30\mu A`$), the shot noise has a slope $`2e`$ as depicted by the dashed line in figure 2. More precisely, the dashed line is : $`S=2eI+\mathrm{1.5\hspace{0.17em}10}^{23}A^2/Hz=\frac{2}{R}(eV+E_c)`$ with $`E_c=22\mu eV`$ and where $`R=V/I`$ is the resistance of the contact. At higher temperature, the general behavior is rounded by thermal fluctuations. Our noise results show two unexpected features: first the crossover to $`S=2eI`$ is predicted for voltage near the gap in the zero temperature models. As said before, in real superconducting-semiconducting junctions exhibiting the reflectionless tunneling a large quasiparticles contribution exists much below the gap. This component explains the observed crossover in accordance with our estimations from conductance measurements, which indicate an increase of the quasiparticles contribution with voltage due to the increase of $`T_{eff}`$ in the silicon layer which is a bad thermal reservoir Quirion00 . The increase of $`T_{eff}`$ could also explain the absence of thermal crossover near $`V5\mu V`$ at $`T=50mK`$, our second unexpected result. However the voltage dependent Johnson noise is rapidly overpassed by the shot noise which grows linearly with the current. The comparison with the noise characteristics at $`T=100mK`$ shows actually that $`T_{eff}`$, which depends in a complicated way on electrical power dissipation via electrons, phonons and superconducting contacts, increases moderately.

In one dimension, Beenakker et al. Beenakker97  calculated the zero energy conductance and shot noise in S-I-N diffusive junctions for various values of $`\mathrm{\Gamma }L/l`$ which corresponds to the ratio of the resistance of the normal metal to the barrier resistance ($`\mathrm{\Gamma }`$ is the barrier transparency, $`L`$ the length of the normal metal and $`l`$ its elastic mean free path). It is shown that, when $`\mathrm{\Gamma }1`$, the shot noise reaches twice the full shot noise for $`\mathrm{\Gamma }L/l0.1`$. At finite energy, Hekking et al. Hekking93  have shown that the length $`L`$ should be replaced by the phase coherence length $`L_\mathrm{\Phi }`$. From ref Quirion00  we get $`0.25`$ for the ratio $`\mathrm{\Gamma }L_\mathrm{\Phi }/l`$ with $`\mathrm{\Gamma }=\mathrm{3.5\hspace{0.17em}10}^2`$, $`L_\mathrm{\Phi }=50nm`$ and $`l=7nm`$. This estimate is in good enough agreement with Beenakker’s predictions since we do not know precisely the properties of the silicon underneath the contactQuirion00 .

In conclusion, we measured the shot noise in a S-I-N junction where a strongly disordered metal (highly degenerate silicon) is in contact with a superconductor (TiN). We found that the shot noise is equal to twice the full shot noise ($`S=4eI`$) at low energy in a regime where electron-hole coherence enhances the conductance. This result is expected from the theory Beenakker97  and corresponds to the Walter Schottky experiment Schottky  with field emission of Cooper pairs through a dielectrics. Electron-hole coherence is required to restore a large enough $`2e`$-(Andreev) component of the subgap current allowing the measurement of the doubled charge. Above $`20\mu V`$, the shot noise follows the Poisson noise $`S_{Poisson}=2eI`$ due to a dominant quasiparticle contribution.

We thank P. Samuelsson, M. Büttiker and X. Jehl for stimulating discussions and J.L. Thomassin for technical assistance.


format: text

text:
# Evolutionary prisoner’s dilemma games with optional participation

## I Introduction

Evolutionary prisoner’s dilemma games (PDG) Hofbauer and Sigmund (1998); Maynard Smith (1982); Weibull (1995); Hamilton (1963); Axelrod (1984) were introduced to study the emergence and maintenance of cooperation among selfish individuals in societies where strategies are either inherited or adopted through basic imitation rules. In the original PDG von Neumann and Morgenstern (1944) two players simultaneously decide whether to cooperate or defect. Mutual cooperation (defection) yields the highest (lowest) *collective* payoff which is shared equally. However, still higher *individual* payoffs are achieved by defectors facing cooperators leaving the latter with the lowest possible payoff. Thus, defection is dominant because they fare better (or at least equal) regardless of the co-players decision. Consequentially, ’rational’ players always end up with the lowest collective payoff instead of the higher reward for mutual cooperation.

Obviously, this result is at odds with observations in human and animal societies. Indeed, cooperation may emerge under certain circumstances (see e.g. Trivers (1971); Nowak and Sigmund (1998); Axelrod and Hamilton (1981); Hauert et al. (2002a); Nowak and May (1992); Sigmund et al. (2001)). In well-mixed populations, i.e. with random matchings, cooperative strategies succeed provided that interactions between the same individuals are repeated with sufficiently high probability. Computer tournaments Axelrod (1984); Nowak and Sigmund (1993) emphasized the importance of particularly simple strategies such as $`D`$ (always defect), $`C`$ (always cooperate), and $`T`$ (”tit for tat”, cooperate on the first move and then repeat the co-player’s move).

In spatially extended systems with limited local interaction, cooperators may thrive by forming clusters Nowak and May (1992); Hubermann and Glance (1993); Nowak et al. (1994) and thereby reducing exploitation by defectors. For example, on square lattices with $`C`$, $`D`$ strategies and suitable parameter values, clusters of cooperators typically expand along straight boundaries while being invaded by defectors along corners and irregular boundaries. The competition between these two invasion processes maintains persistent co-existence of both strategies. When tuning the model parameters, universal phase transitions between states of homogenous defection, co-existence and mutual cooperation are observed Szabó and Tőke (1998); Chiappin and de Oliveira (1999). Very recently the analysis is extended to diluted lattices Vainstein and Arenzon (2001) and to adaptive networks Zimmermann et al. (2001).

The dynamics of such spatial evolutionary PDG becomes increasingly complex if players can adopt three or more strategies Szabó et al. (2001); Killingback and Doebeli (1998); Brauchli et al. (1999). In particular, the crucial role of the $`T`$ strategy is confirmed Nakamaru et al. (1997); Szabó et al. (2000). For example, in an externally driven system supplying $`C`$, i.e. nourishing defectors, cooperative behavior persists through the cyclic dominance of $`D`$, $`C`$, and $`T`$ which maintains self-organizing domain structures Tainaka (1994); Frachebourg and Krapivsky (1998); Frean and Abraham (2001).

In this work, we consider another PDG with three strategies by allowing for optional participation. Players can adopt the $`C`$ and $`D`$ strategies mentioned above or decline participation relying on some small but fixed source of income. Such risk averse players are termed loners ($`L`$).

Voluntary participation in public goods games Hauert et al. (2002a); Szabó and Hauert (2002) turned out to be an efficient way to prevent and avoid the tragedy of the commons Hardin (1968). This system has an inherent cyclic dominance resulting in periodic oscillations in well-mixed populations Hauert et al. (2002b) and self-organizing polydomain structures on square lattices Hauert (2002). However note that the former case requires group sizes larger than pairs i.e. excluding the PDG. Pairwise interactions invariably lead to homogenous states of loners. In sharp contrast, we demonstrate that cooperative behavior persists on square lattices producing self-organizing patterns whereas varying types of oscillations occur on random regular graphs (RRG) emphasizing the importance of topological characteristics. RRGs represent a natural choice for social networks relevant in human behavior as well as economic interactions where the geometrical location hardly matters.

## II The models

We consider the evolutionary PDG in large populations with $`N`$ players in the limit $`N\mathrm{}`$ for different geometries: First, we briefly review the results for well-mixed populations i.e. the mean-field approximation. Second, the players are arranged on a square lattice and interact only with their four nearest neighbors to the north, east, south and west. Finally, the players are located on a RRG with fixed connectivity ($`z=4`$) excluding double or self-connectivities. Therefore, each player has the same number of neighbors as its counterpart on the square lattice. Note that locally a RRG is similar to a tree (or Bethe-lattice) because the average loop size increases with $`N`$ Bollobás (1995). Since boundary problems render Bethe-lattices unsuitable for Monte Carlo (MC) simulations, RRGs serve as suitable substitutes. In addition, the analysis of branching annihilating random walks has justified the validity of pair approximation on RRG Szabó (2000).

Each player follows one of three strategies, namely defection $`D`$, cooperation $`C`$, and loner $`L`$, denoted by a three-component unit vector:

$$𝐬(x)=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right),\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right),\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right).$$
(1)
The total payoff $`m(𝐱)`$ for the player located on site $`𝐱`$ is specified by

$$m(𝐱)=\underset{\delta }{}𝐬^T(𝐱)𝐌𝐬(𝐱+\delta ),$$
(2)
where $`𝐬^T(𝐱)`$ denotes the transpose of $`𝐬(𝐱)`$ and the summation runs over the four nearest neighbors ($`\delta `$) as defined by the geometry under consideration. Following Nowak and May (1992), the payoff matrix $`𝐌`$ is given by

$$𝐌=\left(\begin{array}{ccc}0& b& \sigma \\ 0& 1& \sigma \\ \sigma & \sigma & \sigma \end{array}\right),$$
(3)
where $`b`$ ($`1<b<2`$) determines the temptation to defect. Loners and their co-players invariably obtain the fixed payoff $`0<\sigma <1`$, i.e. they perform better than two defectors but are worse off than cooperative pairs.

Randomly chosen players are allowed to reassess and modify their strategy. For example, the player at site $`𝐱`$ adopts the strategy of one of its (randomly chosen) neighbors $`𝐲`$ with a probability

$$W[𝐬(𝐱)𝐬(𝐲)]=\frac{1}{1+\mathrm{exp}[(m(𝐱)m(𝐲))/K]},$$
(4)
where $`K`$ introduces some noise that occasionally leads to irrational decisions.

The imitation rule or learning mechanism (4) leads to three trivial unstable fixed points, namely the homogenous $`C`$, $`D`$ and $`L`$ states. The highest average payoff ($`m_c=4`$) is achieved if everybody cooperates but this state is vulnerable to exploitation. For sufficiently large $`b`$ the appearance of a single defector suffices for cooperation to break down. The system then evolves into a homogenous $`D`$ state with the lowest average payoff ($`m_d=0`$). Loners may in turn invade this state and thereby provide an escape hatch out of states of mutual defection. Eventually a homogenous $`L`$ state is reached with an average payoff $`m_l=4\sigma `$. Note that solitary defectors are eliminated by noise. Finally, the loop is closed by noting that if no one participates in the PDG, cooperators may thrive again. However, also note that this requires the presence of at least two cooperators. This cyclic rock-scissors-paper type dominance determines the system’s dynamics over a wide parameter range.

## III Mean-field and pair approximations

In the classical mean-field approximation, i.e. in well-mixed populations, this system is described by the (time-dependent) frequencies $`f_i(t)`$ of the three strategies with $`f_C(t)+f_D(t)+f_L(t)=1`$ Szabó et al. (2001); Hauert et al. (2002b). Figure 1 displays the time evolution on a ternary phase diagram.

Regardless of the initial configuration, the system always approaches the homogenous $`L`$ state. This result is robust against changes in the parameter values. The heteroclinic cycle along the boundary of the simplex $`S_3`$ reflects the cyclic dominance of the strategies. Mean-field approximations neglect spatial structures, i.e. they describe individuals with an infinite range of interactions or situations where co-players are randomly chosen.

In contrast, pair approximation is capable of capturing some relevant effects occurring for lasting (fixed) partnership. The frequency of strategies is determined by two-site or pair configuration probabilities (for an introduction see e.g. Marro and Dickman (1999)). The probability of a strategy change depends on the configuration of an eight-site cluster consisting of the randomly chosen player, one of its randomly chosen neighbors (the pair) and their six nearest neighbors. Configuration probabilities of larger clusters are assumed to be products of the corresponding pair configuration probabilities. Due to the same connectivity $`z=4`$ of the square lattice and the chosen RRG, the two structures are indistinguishable by pair approximation.

The following systematic analysis is restricted to the case of $`\sigma =0.3`$ and $`K=0.1`$. Figure 2 shows four typical trajectories for different values of the temptation to defect $`b`$.

For very small $`b`$ ($`1<b<b_0^{(\text{p})}=1.0485(1)`$) the clustering advantage of cooperators suffices to offset exploitation by defectors. Loners are unable to provide a viable alternative and vanish. The ratio of $`D`$ and $`C`$ in the stationary state depends on the parameters $`b`$, $`\sigma `$, and $`K`$. When increasing $`b`$ ($`b_0^{(\text{p})}<b<b_1^{(\text{p})}=1.4670(3)`$), a stable interior fixed point appears, i.e. all trajectories approach this stationary value where all three strategies co-exist. For $`b>b_1^{(\text{p})}`$ periodic oscillations appear. Both, the amplitude and the period increase with $`b`$. In finite populations with $`b>b_2^{(\text{p})}1.85`$ any strategy may go extinct due to noise and the trajectory ends in one of the absorbing homogenous states. Note that in this case the rigorous numerical analysis becomes difficult because of the extremely low configuration probabilities occurring during the oscillations.

## IV Monte Carlo simulations

MC simulations are performed on square lattices and on RRGs. In both cases the runs are started from random initial states. Frequencies and payoffs of all three strategies are traced over time. After an appropriate relaxation time the average frequencies and payoffs are determined (sampling times varied from $`10^4`$ to $`10^5`$ MC steps per site).

The linear extension of the square lattice with periodic boundary conditions was $`400`$ sites for most MC simulations. However, in the vicinity of transition points the system size was increased to up to $`800`$ sites in order to suppress undesired diverging fluctuations.

The simulations confirm that loners die out for small $`b<b_0^{(\text{sq})}=1.0262(1)`$. In fact, since loners invade the territory of defectors, they can only survive if the frequency of defectors exceeds a threshold value. Note that in the vicinity of the transition point, the pair approximation predicts slightly lower defector frequencies and therefore a higher threshold value.

The frequency of loners vanishes continuously when decreasing $`bb_0^{(\text{sq})}`$. The numerical analysis of MC data suggests that this extinction process belongs to the directed percolation universality class Marro and Dickman (1999); Hinrichsen (2000). Related critical transitions were reported and discussed in Szabó and Tőke (1998); Chiappin and de Oliveira (1999); Szabó and Hauert (2002). In contrast to this non-analytical transition, the pair approximation predicts a mean-field type behavior, i.e. the frequency of loners vanishes linearly with $`bb_0^{(\text{p})}`$.

For $`b>b_0^{(\text{sq})}`$ the frequency of loners (defectors) increases (decreases) monotonously with $`b`$. However note that the oscillatory behavior predicted by pair approximation is absent on square lattices. The time evolution of the lattice shows cyclic invasions for any site but the phases of these local oscillations are not synchronized Hauert (2002). Instead, the cyclic invasion sustains a self-organizing poly-domain structure. Similar patterns are observed for spatial rock-scissors-paper games (see e.g. Szabó et al. (2000); Tainaka (1994); Frean and Abraham (2001)).

To illustrate the striking difference between the results obtained on square lattices and RRG, we compare two trajectories determined by MC simulations for the same model parameters (see Fig. 4). On square lattices the trajectory evolves towards a stable interior fixed point. In contrast, on RRG the trajectory approaches a periodic limit cycle. Note that for such large system sizes, the noise amplitude is comparable to the line width. The oscillatory behavior is very similar to the predictions by pair approximation (see Fig. 2). Moreover, the amplitude again increases with $`b`$. To demonstrate this nice agreement, the average, minimal and maximal frequencies of defectors are depicted in Fig. 5.

The MC data is obtained by averaging over $`10^4`$ MC steps per sites after suitable relaxation times for $`N=10^6`$ sites. In absence of periodic oscillations, the difference between the minimum and maximum values increases very slowly with $`b`$ until for $`b=b_1^{\text{rrg}}=1.485(1)`$ a Hopf-bifurcation occurs. For $`b>b_1`$ the extremal values quickly separate from the average value reflecting the oscillating behavior in close agreement with pair approximation.

Simulations for $`b>b_2^{(\text{rrg})}1.82`$ confirm another prediction of pair approximation: in finite populations, oscillations with increasing amplitudes are eventually stopped by the extinction of either one of the strategies and the system ends up in the respective homogenous absorbing state. This means if cooperators are to die out first, the systems evolves into a state with all loners. Obviously any strategy may disappear first but the probabilities depend on the parameters, the initial configuration and the system size. For example, about 87% (13%) of randomly initialized runs lead to a homogeneous loner (defector) state for $`b=1.8`$ and $`N=10^6`$.

On RRG the results of MC simulations are well described by pair approximation as illustrated in Fig. 6.

For such structures, the extinction of loners is expected to be a mean-field type process Szabó (2000). Rigorous analysis of our MC data confirms that loners indeed vanish linearly with $`bb_0^{(\text{rrg})}`$ where $`b_0^{(\text{rrg})}=1.03858(3)`$.

In the region of coexistence ($`b_0^{(\text{rrg})}<b<b_2^{(\text{rrg})}`$) the frequency of loners (defectors) increases (decreases) monotonously with $`b`$. But note that the cooperator frequency in MC simulations slightly increases in the oscillating region ($`b_1^{(\text{rrg})}<b<b_2^{(\text{rrg})}`$). The dashed lines in Fig. 6 represent the unstable stationary solution of pair approximation meanwhile the symbols are determined by averaging over more than 300 cycles in the MC simulations.

## V Conclusions

The present investigations confirm that voluntary participation in PDG interactions is capable of preventing ”the tragedy of the commons” and efficiently avoids states of mutual defection. Introducing the loner strategy gives rise to a cyclic dominance that can maintain the co-existence of all three strategies. With players arranged on a square lattice, the cyclic invasions maintain a self-organizing three-color pattern. In contrast, different types of oscillations in the frequencies of the strategies are observed if players are located on a RRG. Both scenarios have the same connectivity $`z=4`$, i.e. the players invariably interact with four neighbors. Occasionally players reassess their strategy by sampling their neighborhood and switch strategy with a probability depending on the payoff difference.

From the view point of sociology (or behavioral sciences) the most relevant questions refer to the average individual income of defectors, cooperators, and loners. Figure 7 clearly illustrates that spatial extension promotes cooperation.

The average payoff of cooperators is significantly higher than that of defectors but more importantly, the average population payoff is larger than the loners income $`\sigma `$. This means that the possibility to participate in the PDG has a positive net return for the population. At the same time it remains far below the maximum $`m_c=4`$ leaving room for further strategical improvements.

Comparing the results for square lattices and RRG, we note that the average payoff of cooperators, and even more pronounced the average population payoff, turns out to be significantly lower on RRG. Due to the randomly drawn links between the players, cluster formation becomes more difficult while facilitating exploitation. But recall that in absence of spatial structures, only loners survive and forego the chances of the PDG. Similarly, for a sufficiently high temptation to defect ($`b>b_2`$) the oscillations on RRG usually terminate in the absorbing state of all loners (see above).

In this work we investigate effects of the underlying geometry on the success of cooperative strategies by comparing the game dynamics on square lattices and RRG with equal connectivity. These two structures are indistinguishable at the level of pair approximation. However, MC simulations reveal significant qualitative differences arising from the distinct topologies. At the same time, pair approximation turned out to reproduce a detailed picture of the MC results obtained on RRG.

###### Acknowledgements.

G. S. acknowledges support of the Hungarian National Research Fund under Grant No. T-33098. Ch. H. acknowledges support of the Swiss National Science Foundation.


text:
# Coarsening Dynamics of Granular Heaplets in Tapped Granular Layers

## Abstract

A semi-continuum model is introduced to study the dynamics of the formation of granular heaplets in tapped granular layers. By taking into account the energy dissipation of collisions and screening effects due to avalanches, this model is able to reproduce qualitatively the pattern of these heaplets. Our simulations show that the granular heaplets are characterised by an effective surface tension which depends on the magnitude of the tapping intensity. Also, we observe that there is a coarsening effect in that the average size of the heaplets, $`V`$ grows as the number of taps $`k`$ increases. The growth law at intermediate times can be fitted by a scaling function $`Vk^z`$ but the range of validity of the power law is limited by size effects. The growth exponent $`z`$ appears to diverge as the tapping intensity is increased.

In recent years, there has been increasing interest in the collective dynamics of granular materials. The dissipative nature of granular materials gives rise to properties distinct from those of solids and liquids. Many experimental and theoretical attempts have been made to seek a fundamental understanding of this granular state and especially pattern formation in driven granular layers. The experiments include vertically vibrated systems Umbanhowar96 ; bizon98 , tapped or blown thin films of powders Duran00 ; Duran01 , and electrostatically driven granular layers Aranson02 ; Aranson00 . Although many papers have been published to explain the experiments on tapped granular layers, most of them deal with either compactification of thick granular layers Ben-Naim98 ; Knight95 , or static properties of granular heaplets on a thin granular layer Duran01 . The dynamical aspect of the heaping process in a tapped thin granular layer is still not well understood. This paper is concerned with the coarsening dynamics of granular heaplets in tapped granular layers from a theoretical point of view.

This paper is organised as follows. First we introduce a simple model to study this fascinating system by considering the energy dissipation and screening effects during the tapping process. Despite the simplicity, this model is capable of capturing the essential phenomenology, reproducing various morphologies of the coarsening pattern, and showing the way in which heaplets merge as observed in experiments Duran01 . The model is then studied numerically and the results indicate some of the relationships between tapping intensity and the effective surface tension of granular layers. Finally analysis of the results shows that the average size of the the heaplets $`V`$ grows as a power law with the number of taps $`k`$, $`Vk^z`$ for limited range of $`k`$ values. The exponent $`z`$ appears to diverge as the tapping intensity is increased.

Model – Consider a thin layer of $`N1`$ granular particles spread out over a flat plate. The plate is then tapped repeatedly at a low pace with constant shock intensity. In one complete tapping cycle, there are two different processes each of which requires a different description. In the tapping phase, the granular layer is perturbed by an external shock and granular particles acquire kinetic energy to allow them explore the phase space and hop on to neighbouring sites. The hopping range of the granular particles depends on the amount of kinetic energy received by each of the individual particles. The amount of energy supplied to the system is controlled by the dimensionless acceleration $`\mathrm{\Gamma }=U/g\mathrm{\Delta }t`$, where $`U`$ is the velocity of the plate when it is in motion and $`\mathrm{\Delta }t`$ is the time interval over which it is in motion. Note that $`U`$ is also the vertical velocity gained by an elastic particle sitting on the plate when it is tapped. Our particles are not however elastic. Therefore, the vertical velocity gain is $`\alpha U`$, where $`0\alpha 1`$ is the departure coefficient analogous to the coefficient of restitution and it characterises the degree of dissipation of the system. After the tap, the particles undergo a ballistic flight and fall back again onto the static plate where they relax subject to avalanche dynamics before the next tap. In this relaxation phase, the granular particles may stay immobile or move about depending on the local slope. If the local slope is less than a critical slope then the profile remains stationary. Conversely, if the local slope is greater than the critical slope, matter moves down the slope collectively as an avalanche, until the slope is less than the critical slope.

Let the area density $`n(𝐱,k)`$ be the number of granular particles above a unit area of the plate after $`k`$ taps. Here $`𝐱=(x,y)`$, where $`x`$ and $`y`$ are orthogonal coordinates parallel to the plate. The average diameter, $`D`$ of the granular particles is chosen to be $`1`$, so that $`n(𝐱,k)`$ now is equal to the local thickness of the granular layer, so long as there is no compactification throughout the tapping process. Of course this is just an idealisation, the packing fraction of the granular material will be changed Nowak98  if the granular layer is thick or it is intentionally prepared in a low-density state, but such cases will not be considered here.

We can approximate the acceleration $`A_p`$ of the plate’s movement due to tapping as a sum of $`M`$ delta functions

$$A_p(t)=\underset{k=1}{\overset{M}{}}U\delta (t(k1)\tau ),$$
(1)
where the velocity $`U`$ gives the intensity of the shocks which are assumed uniform over the whole plate. $`\tau `$ specifies the time interval between taps. Here $`\tau `$ must be greater than any time scale in the relaxation process, this is to ensure that next tap only occurs after the system is fully relaxed. In the dilute limit, a single inelastic grain sitting on the plate experiences a net force

$$F_z(t)=\alpha \left(\underset{k=1}{\overset{M}{}}mU\delta (t(k1)\tau )\right)mg,$$
(2)
where $`\alpha `$ is the departure coefficient mentioned earlier. During a small interval $`[t_k\mathrm{\Delta }t/2,t_k+\mathrm{\Delta }t/2]`$, where $`t_k`$ is the short hand for $`t_k=(k1)\tau `$, the time of $`k`$-th tap, and the maximum velocity transferred to the particle is

$$\mathrm{\Delta }v_z=_{t_k\frac{\mathrm{\Delta }t}{2}}^{t_k+\frac{\mathrm{\Delta }t}{2}}\frac{F_z(t)}{m}dt=\alpha Ug\mathrm{\Delta }t.$$
(3)
In order for the particle to hop, $`\mathrm{\Delta }v_z`$ must be positive, or $`\alpha U>v_0g\mathrm{\Delta }t`$ so that external acceleration overcomes the gravitational force. In what follows this is always the case.

In the continuum limit one has to make two crucial modifications. $`\alpha `$ must be replaced by an effective departure coefficient of the granular bulk. $`\alpha `$ is expected to be a monotonic decreasing function of $`n(𝐱,k)`$. This is because when the number density is high, inter-grain collisions occur more often and more energy is dissipated, so that grains depart with a smaller departure velocity. While the precise form of $`\alpha (n)`$ is unknown, in this paper it is taken to be

$$\alpha (n)=\{\begin{array}{cc}\frac{\alpha _0\overline{n}}{n(𝐱,k)},\hfill & n>1\hfill \\ \alpha _0\overline{n},\hfill & n1,\hfill \end{array}$$
(4)
for reasons of simplicity and to make comparison with Duran’s model Duran01 . Here, $`\overline{n}`$ is the average density of the system and $`\alpha _0`$ is the effective departure coefficient for average density $`\overline{n}`$. Also, the second term in Eq.(2) needs to be changed. In the continuum limit, according to Duran Duran00  there is an effect dependent on the position of a grain in the heap. If the grain is not at the top of the heap it supports a fraction of the weight of the grains above it, hence its effective mass is increased. As a result, in order for a granular particle sitting on the inclined side of the granular layer to hop, it requires a larger velocity kick than a particle sitting on the flat region of the layer. This screening effect can be incorporated into Eq.(2) by replacing $`m`$ in the second term with an effective mass

$$m^{}=\left((n_Tn)p\mathrm{sin}\theta _c/D\right)m.$$
(5)
Here $`D`$ is the average grain diameter and it is set to 1 henceforth. $`n_T`$ is the altitude of the nearest peak in the corrugated granular layer. $`p`$ is a parameter of unknown value which is set equal to 5 in reference Duran00  and this value appears to give a match to experimental results. $`\theta _c`$ is the angle of repose and is close to the value of $`\pi /6`$Duran00 . After these modifications, Eq.(3) now can be written as

$$\frac{\mathrm{\Delta }v_z}{v_0}=\{\begin{array}{cc}\frac{𝒜}{n}(n_Tn)p\mathrm{sin}\theta _c\hfill & n>1\hfill \\ 𝒜(n_Tn)p\mathrm{sin}\theta _c\hfill & n1,\hfill \end{array}$$
(6)
where $`𝒜=\frac{\alpha _0\overline{n}U}{v_0}`$ is the tapping intensity.

There are close similarities between our model and Duran’s Duran01 . In both models, the resulting pattern formation is due to the competition between the upward hopping motion of particles and the downward screening effect of avalanches. However, the mechanism causing upward motion is different in the two models. Duran Duran01  and Shinbrot Shinbrot98  conjectured that the hopping of granular particles is due to the upcoming air-flux through the porous bed of the granular layer. The velocity of the air-flux at height $`n`$ can be approximated by Darcy’s law $`v_{air}=Kp/n`$, where $`K`$ is the permeability of the granular bed, $`p`$ is the pressure different across the granular bed. The velocity of the air-flux $`v_{air}`$ must be greater than the free fall velocity of a granular particle $`v_f`$ in order to eject this particle. However, in our model we adopt a simpler picture where kinetic energy is transferred to the granular particles via direct collisions between plate and particle and between particle and particle. Energy is dissipated via the effective departure coefficient. The form of the effective departure coefficient is chosen as in Eq.(4) so that the velocity transfer from plate to the particle has the same $`1/n`$ dependency ($`n1`$) as $`v_{air}`$ in Darcy’s law. Obviously, we have the freedom to choose other functional forms of $`\alpha (n)`$, the simple choice here is just to achieve comparison with Duran’s model. As observed from simulation, however, it turned out that the exact form of $`\alpha (n)`$ is unimportant so long as it is a monotonic decreasing function.

Now it remains to specify the hopping process in the tapping phase. When the granular layer is tapped, at an unstable site where $`\mathrm{\Delta }v_z>0`$, particles can hop to any neighbouring site within a circle defined by a maximum horizontal hopping range, $`R_{max}(𝚫v)`$. Here we assume on the unstable site, all the granular particles are re-distributed. Each grain hops with equal probabilities in any direction and with a random hopping range, $`R<R_{max}`$. We can estimate the maximum hopping range $`R_{max}(𝚫v)`$ by the following simple arguments. The maximum velocity received by a particle during one tap is $`\mathrm{\Delta }v_z`$ vertical to the plate. In this case $`R_{max}(𝚫v)`$ is zero for the motion is strictly vertical. It is unlikely that particles will always hop in the vertical direction, an inter-particle collision can change the hopping direction to any angle. Assume that particle is ejected with a speed $`\mathrm{\Delta }v_z\mathrm{cos}\varphi `$ with $`\varphi `$ denoting the angle from the vertical axis of the plate. Then, the horizontal hopping range is

$$R=\mathrm{\Delta }v_z\mathrm{cos}\varphi \mathrm{sin}\varphi t^{}=\frac{2\mathrm{cos}^3\varphi \mathrm{sin}\varphi }{g}(\mathrm{\Delta }v_z)^2,$$
(7)
where $`t^{}`$ is the flight time of the particle in vacuum. $`R`$ is maximum when $`\varphi =\pi /6`$, so that the estimated value for $`R_{max}`$ is given by

$$R_{max}(\mathrm{\Delta }v_z)=\frac{3\sqrt{3}}{8g}(\mathrm{\Delta }v_z)^2=\xi (\frac{\mathrm{\Delta }v_z}{v_0})^2$$
(8)
where

$$\xi =\frac{3\sqrt{3}}{8}v_0\mathrm{\Delta }t.$$
(9)
The estimated $`R_{max}`$ is very crude and does not take into account the details of the air-particle and particle-particle interactions except phenomenologically. However, one can consider $`\xi `$ in $`R_{max}`$ as an adjustable parameter, which varies from one experiment to another experiment depending on the roughness of the material used in experiment and humidity of the surrounding environment.

There are several ways to describe the avalanche process during the relaxation phase. Here we use a slope dependent diffusion equation. We start from the equation of continuity $`_tn=𝐉`$ with the constitutive current $`𝐉`$ given by

$$𝐉=\{\begin{array}{cc}\eta (\beta |n|^21)n,\hfill & |n|>\frac{1}{\sqrt{\beta }}\hfill \\ 0,\hfill & |n|\frac{1}{\sqrt{\beta }},\hfill \end{array}$$
(10)
where $`1/\sqrt{\beta }=\mathrm{tan}\theta _c`$ is the critical slope and $`\eta `$ sets the diffusion rate. The current is chosen in such a way that when the gradient is greater than the critical slope, mass current flows down-hill to smooth out density fluctuations, but no mass is transferred if the gradient is less than the critical slope.

Simulation – The model is studied on a $`N\times N`$ square lattice with a periodic boundary condition. Initially, the granular layer is prepared by assigning a random height between $`[0,5]`$ to each site and letting the granular layer relax. At the beginning, the fluctuation in $`n(𝐱,0)`$ is considerably smaller then the height of the heaplets formed later. On each tapping cycle, the velocity of particles at each site is calculated and the density number is updated according to rules defined by Eq.(6) and Eq.(8). After which the equation of continuity and the corresponding constitutive current equation Eq.(10) are solved and iterated numerically until there is no more mass is transferred at each site.

The following density plots show typical results from the simulation. Fig.1 shows the density plots of $`n(𝐱,k)`$ at different values of $`k`$, the number of taps, for two set of parameters. The top panel shows extended patterns of ridges and corresponds to a smaller tapping intensity ($`𝒜=6.0`$), while the lower panel shows localised heaplets and corresponds to a greater tapping intensity ($`𝒜=16.8`$). Both simulations coarsen when $`k`$ increases, eventually reaching their final stage where the entire system is a single granular heap.

One can see from Fig.1 that as the tapping intensity is increased, the patterns favour isolated circular heaplets which suggests that an effective surface tension can be defined which increases with the tapping intensity. Since there are no attractive forces between granular particles, the granular layer cannot have a true surface tension. This “surface tension” effect is entirely due to the system trying to maximise the local energy dissipation by decreasing the value of the effective departure coefficient so that particles are less mobile. In order to maximise the local energy losses, particles have to be as close to each other as possible, for this will increase the number of inelastic collisions between particles. Therefore, with the constraint that the local slope should not exceeds the critical slope, the global attractor of the pattern is a single large heap, where local energy loss is the largest. Akiyama et al.Akiyama98 , Clément clement92  and Duran Duran01  have discussed the effects of convection on heaping process of granular layers. Duran Duran01  in particular suggests that the “surface tension” is brought about by the convective forces dragging particles from the surroundings into the granular heaplets. However, in our simulation it is due to the maximisation of energy losses, since in our model the granular layer does not contain any information about the internal interactions of the sand heap, therefore convective forces are not taken into account.

Nevertheless, the simulation result does show the Laplace-Young pressure effect suggested by reference Duran01 . A typical example is shown in Fig.2. (The parameters here are $`𝒜=16.0`$ and $`\xi =0.2207`$, but similar results are found over a wide range of values of these parameters.) In the center of the figure there are two connected granular heaplets. After several taps the smaller heaplet merges with the larger heaplet due to matter moving along the connecting neck.

We are also interested in the dynamics of the size of the heaplets. The relevant length scale measure $`l`$ is the average thickness of the site $`h_i`$ weighted by the local volume, $`V_i`$, and the average volume of the heaplet is $`V=l^3`$. Since we know that $`h_i=n_i`$ and the local volume of the site is just $`V_i=n_i\sigma `$, where $`\sigma `$ is the area element of the lattice site, it follows that

$$l=\frac{_ih_iV_i}{_iV_i}=\frac{_in_i^2}{_in_i}.$$
(11)
There is some limited evidence for power law behaviour $`Vk^z`$ in the intermediate region of the log-log plots of $`V`$ versus $`k`$ in Fig.3. There are three regions in each plot. The early region is where the granular layer is randomly distributed, the number density fluctuation is not large enough to trigger coalescence. The late-time region is where all the heaplets have joined into one big heap, leaving some remaining individual particles hopping randomly about and occasionally encountering the large single heap. The intermediate region is also the fast-growing region where heaplets grow and merge into each other. The range of the intermediate region grows and remains linear as the system size increases as can be seen from the two curves with $`𝒜=17.6`$. This shows that system size limits the range of the intermediate power law region. For very large system sizes we would expect on this basis that the intermediate power law region would have a greater range and remain linear in the log-log plot. One can easily estimate the increase in the average volume of the heaplets $`V`$ by the following simple argument. As the saturated value of $`V`$ corresponds to a state where a single heap contains approximately all the grains in the system, it should be proportional to $`\overline{n}L^2`$, here the average number of grains per site $`\overline{n}`$ is a conserved quantity. Then, increment in $`V`$ is given by $`\mathrm{\Delta }\mathrm{log}_{10}V=2\mathrm{\Delta }\mathrm{log}_{10}L`$. In Fig.3 the dashed line corresponds to a lattice size $`L_1\times L_1=80\times 80`$ and the marked-dashed line to $`L_2\times L_2=120\times 120`$, and the vertical shift in the graphs is expected to be $`\mathrm{\Delta }\mathrm{log}_{10}V=0.352`$ which may be compared to $`0.360`$ from the simulation.

The exponent $`z`$ in the power law depends on the tapping intensity $`𝒜`$ as is shown in Fig.4 As we can see $`z`$ is roughly constant, $`z1.0`$, for low tapping intensity. When the tapping intensity is increased $`𝒜>14`$, $`z`$ suddenly increases, and appears to diverge near $`𝒜=18`$. For these large values of the tapping intensity, we observed that the system appears to be in a flat configuration for some time and suddenly a single heap formed within a few taps. For still greater tapping intensities, the heaplet patterns disappear and the layer remains flat with a small random variation superimposed. This pattern-disorder transition occurs at a not very precisely defined critical tapping intensity of $`𝒜18`$. (This critical intensity appears to depend slightly on the initial conditions which is the reason for our statement that it is not precisely defined.) Above this critical tapping intensity, no patterns are observed. This is due to the fact that the perturbation is so strong that no site is stable, and each site constantly undergoes re-distribution on each tapping cycle. Although the power law parametrisation is based on limited evidence, it does provide a useful parameter $`z`$ to distinguish the pattern-forming region ($`𝒜<18`$) from the non-pattern-forming region ($`𝒜>18`$).

We have introduced a simple model of the heaping process in a tapped granular layer. The model is capable of reproducing the essential morphologies of tapped granular systems. Qualitatively the effective surface tension of the granular heaps is closely related to the tapping intensity, and it is shown that there is no need for convective forces for this effective surface tension to exist. The length scale of the system coarsens according to the power law $`lk^z`$ in a limited range. The exponent $`z`$ appears to diverge as tapping intensity is increased and provides a useful parameter for distinguishing the pattern-forming and the non-pattern-forming regions.

We thank Philip Cheung for helpful discussions.


text:
# Spatial Structure of the Internet Traffic

## I Introduction

Internet connects different routers and servers using different operating systems and transport protocols. This intrinsic heterogeneity of the network added to the unpredictability of human practices  make the Internet inherently unreliable and its traffic complex .

There has been recently major advances in our understanding of the generic aspects of the Internet  and web  structure and development. Concerning data transport, most of the studies focus on properties at short time scales or at the level of individual connections , while studies on statistical flow properties at a large scale  concentrate essentially on the phase transition to a congested regime. Despite of these results, large scale studies of traffic variations in time and space are still needed before understanding the new social practices of Internet users.

In this paper, we study the spatial structure of the large scale flow. We present in part II the data studied and in parts III and IV the results of our analysis, showing the existence of a spanning network concentrating the major part of the traffic. Finally, in part V we relate the flow properties and its spatial distribution to scientific activity measured by the number of published papers.

## II Data studied

An important difficulty is to obtain real data measurements of the Internet traffic on a global scale. The availability of data of the French network ‘Renater’ allows us to consider the cartography of Internet’s traffic and its relation with regional socio-economical factors.

The French network ‘Renater’ has about $`2`$ million users and is constituted of a nation-wide infrastructure and of international links . Most of the research, technological, educational or cultural institutions are connected to Renater (Fig. 1). This network enables them to communicate between each other, to get access to public or private world-wide research institutes and to be connected to the global Internet.

We first restrict our analysis to the national traffic and exclude the information exchange with external hosts and routers such as US and Europe Internet or peering with other ISPs. This restriction to a small part of the Renater traffic ($`5\%`$ of roughly $`2000`$ Gigabytes a day) has two methodological advantages: First, it ensures that the traffic studied is strictly professional (mail to non-academics, like family, friends, consultation of newspapers, e-commerce, etc. goes through outer ISP and is not taken into account); Second, it helps to understand the regional traffic structure and its relation with local economical factors. We believe that the global patterns emerging for the Renater network will be relevant for larger structures such as the global Internet.

The data consist of the real exchange flow (sum of Ftp, Telnet, Mail, Web browsing, etc.) between all routers, even if there is not a direct (physical) link between all of them. For a connection $`(i,j)`$ between routers $`i`$ and $`j`$ ($`ij`$), $`F_{ij}(t)`$ (in bytes per $`5`$ minutes) is the effective information flow at time $`t`$ going out from $`i`$ to $`j`$. For technical reasons, data for a few routers were not reliable and we analyzed data for $`26`$ routers which amounts in $`26\times 25`$ matrices $`F_{ij}(t)`$ given for every $`\mathrm{\Delta }t=5`$ minutes for a two weeks period (the quantities $`F_{ii}`$ are excluded from the present study).

As an example of the measured time-series, we show (Fig. 2) the information flow versus time between two routers located in Grenoble and Marseille for a nine days period. One can see the different days and within days, bursts of intense activity. In this study, we focus on the flow and not on the growth rate (and its correlations) used in a previous study .

## III Databases versus Active Centers

We now present our empirical results. The time averaged incoming flow

$$F_{in}(i)=\underset{j}{}\overline{F}_{ji}$$
(1)
at a given router $`i`$ is a measure of the Internet activity of the corresponding region (the over-bar denotes the average over time). On the other hand, the average outgoing flow

$$F_{out}(i)=\underset{j}{}\overline{F}_{ij}$$
(2)
can be interpreted as the total request emanating from other routers. It is thus a measure of the degree of interest produced by this region.

We plotted both quantities $`F_{in}`$, $`F_{out}`$ versus their rank (Fig. 3a,b). In contrast with many cases observed since the work of Zipf , the observed distributions are not power laws but exponentials. This might be the signature of a transient regime and would mean that the Internet didn’t reach his stationary state, but it is more probably the sign that the Internet traffic has a unique, non hierarchical-type structure . This exponential behavior also means that–at least in the Renater network–there are essentially two categories of regions. Considering Internet activity (Fig. 3a), one can distinguish active from (almost) inactive regions. Roughly, there are about eight cities which receives $`80\%`$ of the total traffic, the rest being (exponentially) negligible. Concerning the outgoing flow (Fig. 3b), there are about five most visited regions, the rest being comparatively ‘unattractive’. We checked that for different time windows the order of these cities can slightly change, but the exponential behavior is independent of these seasonal effects.

It is interesting to note that the most active and visited regions are not the same showing that each region has its specific activity. Regions with a large incoming flow can be classified as active research centers with a great need of information, and regions with a large outgoing flow correspond to important information resources such as e.g., databases or libraries.

At this stage, we have shown that in the Renater traffic there is a small number of receivers (located in active regions) and emitters (visited databases). However, a further question concerns the secondary routers and the fine structure of flows. Indeed, $`F_{in}`$ (and similarly for $`F_{out}`$) could be a sum of many small contributions coming from many regions, or in contrast there could be only few regions which exchange a significant flow. Simple quantities which can characterize the fine structure of the incoming flow at router $`i`$ are the $`Y_k`$’s introduced in another context

$$Y_k(i)=\underset{j=1}{\overset{N}{}}(W_{ji})^k$$
(3)
where $`W_{ji}=F_{ji}/F_{in}`$ is the weight associated with the incoming flow $`F_{ji}`$ (and similar expressions for the structure of outgoing flow). It is easy to see that $`Y_0=N`$, $`Y_1=1`$ and the first non trivial quantity is $`Y_2`$. We can illustrate the physical meaning of $`Y_2`$ with simple examples. If all weights are of the same order $`W_{ji}1/N`$ for all $`i,j`$ then $`Y_21/N`$ is very small. In contrast, if one weight is important for example of the order $`1/2`$ and the others negligible $`1/2(N1)`$ then $`Y_21/4`$ is of order unity. Thus $`1/Y_2`$ is a measure of the number of important weights. We plot $`Y_2`$ for both the incoming and outgoing flows (the statistics is over two weeks). The result (Fig. 4) shows clearly that the most probable value is $`1/4`$ and that $`Y_2`$ is larger than $`1/N1/300.03`$ (except for few cases which appear in the histogram). This confirms the fact that a few routers are exchanging most of the information, the rest of the network being negligible.

## IV Spanning network

In order to illustrate the above results, we construct the network $`S_k`$ connecting a number $`k`$ of routers and carrying the maximal flow denoted by $`F(S_k)`$. We increase $`k`$ from $`2`$ to $`26`$ and we obtain the result plotted in Fig. (5a). It appears clearly that a small fraction of links carry most of the flow. This behavior is encoded in the fact that $`F_{tot}`$ is a power law for $`p0`$

$$F_{tot}(S_k)p^\theta $$
(4)
with an exponent smaller than one ($`\theta 0.6`$), so that a small variation of the number of connections leads to a large variation of the transported flow.

This analysis completes the Renater traffic map: there is a small number of receivers and emitters exchanging significant information between them, the rest of the network being exponentially negligible. This demonstrates the existence of a ‘spanning network’ carrying most of the traffic and connecting the main emitters to the main receivers.

In order to visualize this spanning network on the French map, we filter the flow with the following procedure. We first select flows above a certain threshold $`F_c`$, and then we select a connection $`(ij)`$ only if the corresponding flow $`F_{ij}`$ represents a large percentage of (i) the outgoing flow from $`i`$ and (ii) the incoming flow in $`j`$. The result is shown on Fig. 5b. We checked that the instantaneous (or averaged over a different time window) spanning network is eventually different, but always with the same characteristics (small number of interconnected emitters and receivers). The procedure described above could thus be used as a simple filter in order to visualize in real-time a complex flow matrices.

## V Digression to Scientometrics: Internet Traffic versus Scientific activity

So far, we have studied statistical properties of the traffic, but an important point is to relate them to economical or social factors. The Internet activity should in principle be related to social pointers such as the number of inhabitants, the number of students, and so on. From our data, the indicator which shows the best correlation with Internet activity is scientific activity measured by the number of published papers . One can expect that the more a scientist consults books or data, the more he/she will publish. This principle, ‘the more you read, the more you write’, although commonly accepted in a number of historical cases , is difficult to evaluate quantitatively. The main difficulty being the measure of the amount of information gathered by scientists in libraries. In the case of Internet, the information gathered by scientists working in a given region can be estimated by the average incoming flow in the corresponding router. Since the information needed for a scientist is usually scattered world-wide, it is important here to take into account the total incoming flow, including exchanges with international hosts. We thus compare the total average incoming flow (per scientist) with the average number of papers published (per scientist) per year by the region’s universities (obtained from the SCI database). As a representative panel, we choose to use data about papers published only by scientists in the national research institution (CNRS ). We represent these data on Fig. 6. This plot shows that the average incoming flow per scientist $`\overline{f}`$ in a region is increasing with the number of published scientific papers per scientist $`p`$ by this region’s laboratories roughly as a power law

$$\overline{f}p^\beta $$
(5)
with exponent $`\beta 1.1\pm 0.1`$. This result confirms quantitatively the intuitive principle stated above and is particularly interesting from the point of view of the Web’s social impact. Indeed, it implies that the number of publications is growing with the incoming flow as a power law with exponent $`1/\beta 1`$: the more one uses Internet, the more one publishes! This result indicates that on average the use of Internet has a positive impact on research productivity.

## VI Conclusion

In summary, we have shown that the major part of the traffic takes place only between a few routers while the rest of the network is almost negligible. We have proposed a simple procedure to extract this (bipartite) spanning network, which could have some implications in visualization and monitoring of real-time traffic. In addition, resources allocation and capacity planning tasks could benefit from the knowledge of such a spanning network. These results point towards new ways of understanding and describing real-world traffic. In particular, any microscopic model should recover these statistical properties and our results provide a quantitative basis for modeling the dynamics of information flow.

We also have shown that the scientific activity of a region is increasing with its Internet activity. This indicates that it is difficult for a scientist to avoid the use of Internet without affecting his/her productivity measured in terms of publications. This result also demonstrates that in addition to increase people’s social capital  the Internet has a measurable positive impact on research production. More generally, it underlines the importance of Internet as knowledge sharing vector. This study also suggests that the Internet activity could be used as an interesting new socio-economical pointer well adapted to the information society.

Finally, these results exhibit some global statistical patterns shedding light on the relations between the Internet and economical factors. It shows that in addition to the structural complexity of the web and the Internet, the traffic has its own complexity with its own cartography.


text:
# Untitled Document

Write to:
adelpop@unibg.it in order to get an electronic copy of the paper.


text:
# Selective advantage of topological disorder in biological evolution

## 1 Introduction

The mechanism of biological evolution is a very challenging topic for the physical community. This is well expressed in the numerous models of biological evolution that have emerged in recent years. The list of models studied starts with large-scale properties of the evolutionary process—massive global extinctions—and ends with the works which are aimed at following the replicative behavior of the chemical structures holding information, namely DNA molecules. A lot of effort has been devoted to this area recently ba\_fly\_la\_92 ; ba\_sne\_93 ; pa\_ma\_ba\_96 ; va\_au\_95 ; fe\_pla\_dia\_95 ; so\_ma\_96 ; ro\_ne\_96 ; zhang\_97 ; drossel\_98 ; sol99 ; sla\_kot\_99 ; sla\_kot\_00 ; kot\_sla\_sta\_02 ; nim\_cru\_huy\_99 ; wilke\_01 ; wilke\_01a  and still new fruitful ideas emerge kra01 ; chri\_coll\_hall\_jen02 ; hall\_chri\_coll\_jen02 .

The process of biological evolution consists in three steps: $`reproductionmutationselection.`$ The important thing to note is that the *biological fitness* function which denotes an individuals’ ability to produce viable offsprings depends on their phenotype. On the other hand, the mutations occur in the genotype—information stored in a sequence of DNA. The fitness function, which assigns to each microscopic genotype its ability to reproduce and pass on to the next generation, shows an overwhelming complexity. This property makes the theoretical treatment of evolution an extremely complicated task.

Simplifications of the problem are necessary. The fruitful idea of adaptive landscapes was introduced by Wright wri31  and was later further simplified to fitness landscapes later. In this scheme the individual is represented as a point in a multidimensional space and to every point a fitness value is assigned. Thus, a landscape is formed of mountains of genetically adapted positions and valleys of lethal genomes. Note the idealization: the fitness is directly given by the genome of the individual, not by its (extended) phenotype.

In fact, the fitness landscape is not static since it strongly depends on the ever changing environment, which includes interactions with (co)evolving species as well as abiotic influences kauffman\_90a ; kau91 ; wil01 . Evolutionary process manifests itself the ascent of individuals to peaks on the fitness landscape.

A broad set of different fitness landscapes was recently used to study the behavior of the evolutionary system. These models employ both static pel01 ; alt01  and dynamic landscapes wil01 ; nil00 ; pekalski\_00 ; pek\_wer\_01 . They include the *sharply-peaked landscape* (SPL) with a single preferred genome (wildtype), the *Fujiyama landscape*, or the *holey landscape* (HL). Several recent reviews summarize various approaches explored wil01 ; baa99 ; pel97 ; drossel\_01 .

Generally, one simplifies the genetic code considering only a two-letter $`\{0,1\}`$ alphabet. Let us consider, for example, the digits $`1`$ and $`0`$ as symbols for two different alleles of a certain gene111Other possibility is to assign pyrimidines in the genetic code by $`1`$ and purines by $`0`$.. Assuming constant genome length of $`d`$ loci the state space is a hypercube $`𝒢=\{0,1\}^d`$.

The dynamics of evolution on fitness landscapes have the most prominent scheme in the *quasispecies* model as introduced by Eigen eig71 . Originally it was introduced as a model for chemical prebiotic evolution, but it showed itself plausible in the investigation of the mechanisms of microevolution of viruses and bacteria, i.e. organisms with relatively simple genomes.

The quasispecies is defined as a cloud of closely related individuals. They hang together but certainly they may move on the fitness landscape. We assume an infinite population size and therefore the dynamics of the quasispecies are probabilistic, but without any noise which would be induced by finite size effects. The quasispecies obeys three basic processes: reproduction, selection and is liable to mutational genetic changes. Reproduction and selection are treated together—the reproductive ability depends on the fitness and hence selection takes part here. Mutation can occur in the individual’s genome with a probability rate $`\mu `$.

It is evident that the native geometry of the genome space is the above mentioned hypercube $`𝒢=\{0,1\}^d`$. So, this is the natural starting point when building a model for a fitness landscape. Recently, the main focus has been aimed at the presence or absence of the adaptive regime induced by a single maximum in the fitness landscape. Therefore, the first thing to try is the sharply-peaked landscape on the hypercube. This model was solved exactly by Galluccio et al. ga\_gra\_zha\_96 ; galluccio\_97 .

The most important feature found in the quasispecies model in SPL is the error threshold tar92 ; fra97 . It is the phase transition that separates two regimes of the quasispecies evolution, namely the adaptive regime and the wandering regime. In the adaptive regime the localized cloud of quasispecies is formed around the wildtype (for us a certain site in the lattice), whereas in the wandering regime no quasispecies is formed due to a mutation load.

The error threshold phenomenon comes out of the competition between the selective advantage of a certain wildtype and the mutation load presented by the rate $`\mu `$. The threshold is then characterized by the specific value of the selective advantage. We want to show that the selective advantage in the specific site is not the only parameter, whose value can distinguish two significantly different regimes. There are also geometrical properties, e.g. connectivity of the site, that can make the genotype in the specific site advantageous. These are the main tasks of this article.

Indeed, the real landscapes are far more complicated than SPL. It is expected that rugged landscapes with many competing maxima represent a realistic picture ba\_fly\_la\_92 . Such landscapes are well-known in the theory of spin glasses me\_pa\_vi\_87  and a “spin-glass” theory of evolution was investigated, e.g. in am\_pe\_sa\_91 . In our previous work we have studied several similar models of the fitness landscape, too kol\_sla\_02 .

Yet another approach to the modelling of the fitness landscape is used for computations gav98 ; gav99 —the so-called holey landscape (HL), where the fitness is either a positive constant (which may be set to $`1`$) or it is $`0`$ which means that the individual with the corresponding genetic code dies with a probability of $`1`$ without having offspring.

Indeed, a large part of the point mutations which may occur at the basic level of the evolutionary picture are lethal for the individual—they lead to the ‘lethal’ sites. Therefore, the hypercube does not represent a good approximation to the evolutionary dynamics, because only a small part of its edges represent paths to possible new ‘hospitable’ genomes.

A sparsely connected set of points selected at some of the hypercube corners is perhaps a better choice. Such sparse graphs and their adjacency matrices are under study extensively at present; see e.g. bi\_mo\_99 ; monasson\_99 ; sem\_cug\_02  and references therein. The authors observe the localization of the eigenvectors of the sparse matrices due to topological characteristics only.

We approach the problem from a somewhat different perspective which resembles the study of topologically disordered solids where the random network of bonds is often well modelled by the Bethe lattice wil\_mas\_vel\_86 . This, of course, supposes that there are no short loops in the graph. Supposing we are above the percolation threshold, our random lattice forms a giant cluster within the hypercube and the typical length of loops is in the number of sites $`N=2^d`$ of order $`\mathrm{log}N`$, where $`N`$ is the total number of sites $`N=2^d`$, see bollobas\_85 . Thus, the Bethe lattice may be a good model for the topology of our sparse graph.

The main question addressed in our work will be the following: Selective advantage in biological evolution is usually attributed to higher reproductive success. If advantage due to high individual fitness exceeds a certain threshold, a quasispecies is formed around the site. This is the usual error-threshold phenomenon. Here we ask, whether some factors related purely to the structure or topology of the genome space may lead to similar selective advantage, and if a certain threshold can be found separating the adaptive and wandering regime of biological evolution.

## 2 Model

Recently, kol\_sla\_02  we have modelled evolution in a holey landscape using the regular Bethe lattice. Now we will try to represent it in a more precise way and take into account the irregularity of the lattice. For each site we select whether the links that leave it lead to another ‘hospitable’ site of the hypercube, or not. In the latter case, there are some links leading from the site to the ‘lethal’ sites of the hypercube and so all mutants that took this direction are doomed.

We start with a regular Bethe lattice with a connectivity $`k`$. We construct the irregular Bethe lattice by randomly assigning lethal sites and removing all sites which are connected to the rest of the lattice only through a lethal site. The evolution process amounts to diffusion on this lattice. Hops from a viable site to any of its neighbors occur with an equal rate $`\mu `$. Hops from lethal sites are prohibited. So, the edges to lethal sites are “dead ends” or “dangling bonds”, in the language of condensed matter physics.

Let $`i`$ be a viable site. Let us denote $`(i)`$ as the set of viable neighbors of $`i`$. $`\kappa _i=|(i)|k`$ is the number of those viable neighbors. For $`k`$ we suppose only that it is large enough to fulfill the previous inequality. The arbitrariness of the choice of $`k`$ will be clarified later on. Intuitively, during evolution the probability will flow out from $`i`$ to all its $`k`$ neighbors, but flow in only from $`\kappa _i`$ viable ones. Therefore, we expect that a site with larger $`\kappa _i`$ is more likely to gather individuals and form a quasispecies cloud. The scope of this paper is to elaborate this intuitive picture in a more formal manner.

If we denote $`p_i(t)`$ as a (relative) population222By “population” we mean the infinite population limit, the probability of finding a given individual in the site $`i`$ of the lattice. We prefer the term “relative” population because of the perturbations that will be added to the lattice and will destroy the conservation of the overall relative population. of site $`i`$ at time $`t`$, we can write the following master equation

$$\dot{p}_i(t)=\underset{j}{}T_{ij}p_j(t)$$
(1)
where the matrix $`T`$ contains the effects of mutations and reproduction:

$$T_{ij}=\{\begin{array}{cc}\hfill \mu k+\zeta & i=j\hfill \\ \hfill \mu & j(i)\hfill \\ \hfill 0& \text{elsewhere.}\hfill \end{array}$$
(2)
The constant $`\zeta `$ is introduced artificially in order to keep the total population constant. In fact, there are two causes of the net population outflow which must be countered by the $`\zeta `$ term. The first one is the presence of traps—lethal sites that absorb the individuals. The second one comes from the very topology of the Bethe lattice. Indeed, any finite Bethe lattice has the rather counter-intuitive property that the number of the surface sites is comparable to the number of bulk sites and this property holds even in the limit of an infinite number of sites. Therefore, there is always a net flow of probability toward the surface. However, we are interested in the properties of the sites deep in the bulk and the outflow toward the surface should be considered as an artifact of the Bethe lattice approximation. To sum up, we will fix the value of $`\zeta `$ later on in the calculations, in order that the population remains fixed.

The dynamical matrix $`T`$ is in the model symmetric, and this enables us to use the method of the resolvent in our calculations. The symmetry results from the assumption of equal mutation rates, $`\mu `$. This widely used approximation simplifies all the following calculations and, since we are interested only in a stationary state of the system, it is not considered to change the main results. One can introduce the edge-dependent rates $`\mu _{ij}`$ and then $`\zeta `$ would become a function of all $`p_i(t)`$ and Eq. (1) would generally become non-linear one. But this approach goes beyond the scope of our article.

## 3 Partitioning

As in the previous paper kol\_sla\_02  we investigate the formation of a localized state, now interpreted as a quasispecies cloud, through the properties of the resolvent of the matrix $`T`$

$$𝒢(z)=(zT)^1.$$
(3)
The idea of the calculation is quite simple: In the long-time limit only the largest eigenvalue of the matrix $`T`$ survives, and the corresponding eigenvector describes the stationary state of the evolutionary system.

In order to find the largest eigenvalue of $`T`$, we search for the poles of an element of the resolvent matrix $`𝒢`$. These are exactly the eigenvalues of $`T`$, no matter which element of $`𝒢`$ we choose. And since we want to observe the formation of the localized state around a specific site $`i=1`$, we need the diagonal matrix element of the resolvent

$$G(z)=[𝒢(z)]_{11}$$
(4)
which can be calculated using the partitioning (projector) method lowdin\_62 , explained in more detail in kol\_sla\_02 . The loop-less structure of the Bethe lattice greatly simplifies the treatment. We proceed essentially in two steps, which are illustrated in the Fig. 1. In the first step, we project out of the site $`i=1`$ itself. The rest of the Bethe lattice splits into disconnected branches. We find

$$G(z)=\frac{1}{z+\mu k\zeta \mu ^2_{j(1)}\mathrm{\Gamma }_j(z)},$$
(5)
where $`\mathrm{\Gamma }_j(z)`$ is the diagonal element of the projected resolvent on the terminal site of the branch starting with site $`j`$. (The terminal sites are denoted as $`2`$ in Fig. 1.)

In the second step, we calculate $`\mathrm{\Gamma }_j(z)`$ by projecting out the site $`j`$. We find, similar to the previous case,

$$\mathrm{\Gamma }_j(z)=\frac{1}{z+\mu k\zeta \mu ^2_{l(j)\backslash \{1\}}\mathrm{\Gamma }_l(z)}.$$
(6)
Some of the sites $`l(j)\backslash \{1\}`$ are denoted by $`3`$ in Fig. 1. Iterating the equation (6) we could in principle calculate the resolvents on the terminal sites and inserting them into Eq. (5) obtain the desired quantity.

This procedure works well in an infinite regular Bethe lattice, where $`\mathrm{\Gamma }_j(z)`$ for $`j`$ deep in the bulk does not depend on the site index $`j`$ and (6) represents in fact a closed equation for $`\mathrm{\Gamma }(z)\mathrm{\Gamma }_j(z)j`$. However, this procedure cannot be directly applied in the case of an irregular Bethe lattice. Nevertheless, it is a good starting point for an approximation we will describe in the following.

## 4 Effective medium approximation

The mean-field-type treatment of disordered solids was developed a long time ago within the coherent potential approximation (CPA) vel\_kir\_ehr\_68 . In the theory of sparse random matrices it was elaborated using the replica method bi\_mo\_99 ; monasson\_99 ; sem\_cug\_02  and called the effective medium approximation (EMA). In this section we will use the EMA for the irregular Bethe lattice without using the replica trick.

The main idea relies on a simple observation that the sum $`_{l(j)\backslash \{1\}}\mathrm{\Gamma }_l(z)`$ containing $`\kappa _j1`$ terms can be replaced by its average value for large $`\kappa _j`$, thus neglecting fluctuations. Indeed, it was proved that the CPA is exact in the limit of infinite connectivity jan\_vol\_92 . Therefore, our version of the EMA amounts to approximating

$$\underset{l(j)\backslash \{1\}}{}\mathrm{\Gamma }_l(z)(\kappa _j1)\mathrm{\Gamma }(z).$$
(7)
In order to close the equations, we must average expression (6) over the probability distribution of connectivities $`P(\kappa )`$

$$\mathrm{\Gamma }(z)=\underset{\kappa }{}\frac{P(\kappa )}{z+\mu k\zeta \mu ^2(\kappa 1)\mathrm{\Gamma }(z)}.$$
(8)
The latter equation (8) is the core of the EMA. When we insert its solution $`\mathrm{\Gamma }(z)`$ into (5) and perform again the average over connectivities, we obtain the averaged diagonal element of the resolvent $`G(z)`$, which is now site-independent. The parameter $`\zeta `$ represents the shift in the variable $`z`$ and should be adjusted so that the upper edge of the support of the imaginary part of $`G(z)`$ (i.e. the density of states) lies at $`z=0`$. This expresses the requirement of the conservation of the population size.

As an illustration we show in Fig. 2 the real and imaginary parts of $`\mathrm{\Gamma }(z)`$ for the connectivity distribution chosen as

$$P(\kappa )=p\delta (\kappa 2)+(1p)\delta (\kappa 3)$$
(9)
for $`0p1`$. We can see that the imaginary part of $`\mathrm{\Gamma }(z)`$ approaches zero as $`z^{1/2}`$ at the band edge and the real part approaches a finite limit. From the technical point of view, the finiteness of the limit is the source of the transition between localized and delocalized states.

## 5 Single defect

So far we have investigated the Bethe lattice as an averaged homogeneous effective medium. Now we investigate the behavior of the resolvent at a site $`i=1`$ with a specific connectivity $`\kappa _1`$. The situation is schematically depicted in Fig. 3.

This approach has a biological motivation. What we see in nature is that there are frequent groups of closely related animals. These groups form taxonomic classes or, better said, the classes are defined as the groups of such closely related individuals. In the following, we examine if the clustering of individuals, modelled as the addition of the site with largest connectivity, leads to some observable changes in the quasispecies evolutionary process.

The on-site resolvent corresponding to this site can be found from Eq. (5), where the resolvents at the terminal sites are approximated as in (7). This is the essence of the single-defect approximation (SDA). Hence

$$G_{\mathrm{SDA}}(z)=\frac{1}{z+\mu k\zeta \mu ^2\kappa _1\mathrm{\Gamma }(z)}.$$
(10)
A state localized at this site exists, if the resolvent $`G_{\mathrm{SDA}}(z)`$ has a pole for a real positive $`z`$. This is equivalent to the condition

$$\kappa _1>\kappa _c\frac{\mu k\zeta }{\mu ^2}\frac{1}{\mathrm{Re}\mathrm{\Gamma }(0)}.$$
(11)
Therefore, the state is localized and the quasispecies cloud is formed for integer $`\kappa _1`$ greater than $`\kappa _c`$. This corresponds to the adaptive regime of the evolutionary process. We use as an example again the distribution (9) and show in the Fig. 4 the dependence of the threshold $`\kappa _c`$ on the parameter $`p`$ of the probability distribution.

We will consider yet another type of “defect” in the effective medium. Imagine that the site investigated, $`i=1`$, corresponds to a genome of length different to the average. This difference is translated in our model as a choice of different total connectivity $`k`$ in the site $`i`$ (including edges to both viable and lethal sites). So, we may replace $`kk+\lambda `$ in Eq. (10). Then we may investigate the transition from the wandering to the adaptive regime when varying two parameters $`\lambda `$ and $`\kappa _1`$. We will see that the individuals added due to the $`\zeta `$ term in Eq. (1) are redistributed in the lattice in such a manner that the total population changes. This is only due to addition of the “defect” and the localization of individuals in its vicinity. The condition for the existence of a localized state, i.e. for the adaptive regime, is analogous to (11) and can be written as

$$\lambda <\lambda _c\mu \kappa _1\mathrm{Re}\mathrm{\Gamma }(0)\frac{\mu k\zeta }{\mu }.$$
(12)
We show in Fig. 5 an example of the dependence of $`\lambda _c`$ for the same connectivity distribution (9) as in previous cases, and for a specific, fixed choice of $`\kappa _1`$.

## 6 Conclusions

We modelled the complex fitness landscape of biological evolution with an irregular Bethe lattice. The formation of localized quasispecies in the adaptive regime was observed via the occurrence of the isolated pole in the on-site resolvent. Using the effective medium approximation we calculated the disorder-averaged resolvent and within the single-defect approximation we investigated the localization.

We found that the quasispecies is formed if the site is connected to a sufficiently large number of viable sites. We found the condition for the critical value of the connectivity $`\kappa _c`$ separating the adaptive evolutionary regime for $`\kappa >\kappa _c`$ and the wandering regime for $`\kappa <\kappa _c`$. As expected, $`\kappa _c`$ grows with the average connectivity and qualitatively speaking the quasispecies are formed at such sites that have sufficiently larger number of viable neighbors than average. We may interpret it as a purely topological selective advantage: The species has better chance to survive not because of its individual reproductive abilities, but because it is less vulnerable to random alterations of the genetic code.

We investigated also another deviation from the mean, connected to the overall number of neighbors. In more biological language it corresponds to the length of the genetic code. We observed that the adaptive regime is favored for lower values of the connectivity, i.e. for shorter genomes, if the number of viable neighbour sites is considered constant. This type of selective advantage is therefore also of topological origin and has a similar biological interpretation to that presented above. Indeed, longer genome means higher probability of a lethal mutation.

In both cases it is important to note that it is the deviation from the average topology which makes the selective advantage work and which leads to the formation of localized states. So, it is the (sufficiently strong) topological disorder that is responsible for the formation of quasispecies.

We may conclude by summarizing that without resort to individual reproductive capacities, biological evolution favors genomes which are shorter and more robust to random mutations. This has one more important implication; a genome, which can be easily mutated without affecting the death of its carrier, means also a less-defined species. One may therefore predict that successful species will exist in a broad variety of slightly different sub-species. This effect is caused by the selective advantage of certain topologies of the genome space. Observation of the variability within a single species may therefore say something of the relative importance of topological selective advantage in comparison to individual reproductive success.

###### Acknowledgements.

The authors would like to thank to Jan Mašek for critical reading of the manuscript. Michal Kolář would like to thank to Vladislav Čápek for fruitful discussions and to Anton Markoš for helpful discussions on the topic of biological evolution, too.