Universality in survivor distributions:Characterising the winners of competitive dynamics

Universality in survivor distributions:
Characterising the winners of competitive dynamics

J. M. Luck jean-marc.luck@cea.fr Institut de Physique Théorique, Université Paris-Saclay, CEA and CNRS, 91191 Gif-sur-Yvette, France    A. Mehta anita@bose.res.in S. N. Bose National Centre for Basic Sciences, Block JD, Sector 3, Salt Lake, Calcutta 700098, India

We investigate the survivor distributions of a spatially extended model of competitive dynamics in different geometries. The model consists of a deterministic dynamical system of individual agents at specified nodes, which might or might not survive the predatory dynamics: all stochasticity is brought in by the initial state. Every such initial state leads to a unique and extended pattern of survivors and non-survivors, which is known as an attractor of the dynamics. We show that the number of such attractors grows exponentially with system size, so that their exact characterisation is limited to only very small systems. Given this, we construct an analytical approach based on inhomogeneous mean-field theory to calculate survival probabilities for arbitrary networks. This powerful (albeit approximate) approach shows how universality arises in survivor distributions via a key concept – the dynamical fugacity. Remarkably, in the large-mass limit, the survival probability of a node becomes independent of network geometry, and assumes a simple form which depends only on its mass and degree.

I Introduction

The fate of agents in any situation where death is a possibility attracts enormous interest: here, death does not have to be literal, but could also refer to bankruptcies in a financial context or oblivion in the context of ideas. This panorama of situations is usually captured by agent-based models including predator-prey models, the best known being the Lotka-Volterra model Murray (1989); Berryman (1992); Szabó and Fáth (2007) or indeed, winner-takes-all models, which embody the survival of the largest, or the fittest. The model we study here belongs to the second category; the most massive agents grow at the cost of their smaller neighbors, which eventually disappear. Motivated by the physics of interacting black holes in brane-world cosmology Majumdar (2003); Majumdar et al. (2005), its behavior in mean-field and on different lattices was investigated at length in Luck and Mehta (2005), followed by simulations on more complex geometries Thyagu and Mehta (2011); Aich and Mehta (2014). The essence of the model in its original context Majumdar (2003); Majumdar et al. (2005) involved the competition between black holes of different masses, in the presence of a universal dissipative ‘fluid’. The model also turned out to have a ‘rich-get-richer’ interpretation in the context of economics, where it was related to the survival dynamics of competing traders in a marketplace in the presence of taxation (dissipation) Mehta (2013); Mehta et al. (2005).

One of the most important questions to be asked of such models concerns the distribution of survivors, i.e., those agents who survive the predatory dynamics. In spatially extended models, the presence of multiple interactions makes this a difficult question to answer in full generality. The seemingly universal distributions of survivor patterns put forward in Thyagu and Mehta (2011) motivated us to ask: can one provide a theoretical framework for the appearance of some universal features in survivor distributions? The answer is yes, as we will demonstrate in this paper.

We first introduce a much simpler version of the model investigated earlier Luck and Mehta (2005); Thyagu and Mehta (2011); Mehta (2013); Aich and Mehta (2014), which, however, preserves features such as the exponential multiplicity of attractors (Section II) essential to its complexity. We characterize exactly the attractors reached by the dynamics on chains and rings of increasing sizes (Section III), an exercise which illustrates how a very intrinsic complexity makes such computations rapidly impossible. This leads us to formulate an approximate analytical treatment of the problem on random graphs and networks (Section IV), based on inhomogeneous mean-field theory, which yields rather accurate predictions for the survival probability of a node, given its degree and/or its initial mass.

Ii The model

The degrees of freedom of the present model consist of a time-dependent positive mass at each node of a graph. These masses are subject to the following first-order dynamics:


Each node is symmetrically (i.e., non-directionally) coupled to its neighbors, which are all nodes connected to  by a bond; is a positive coupling constant.

The dynamics are deterministic, so that all stochasticity comes from the distribution of initial masses . This model is a much simpler version of the one investigated in Luck and Mehta (2005) which was inspired by black-hole physics Majumdar (2003); Majumdar et al. (2005). In particular the explicit time dependence and initial big-bang singularity of that model are here dispensed with; only a simple nonlinearity, quadratic in the masses, is retained. Despite these simplifications, the present model keeps the most interesting features of its predecessor, such as those to do with its multiplicity of attractors.

The coupling constant can be scaled out by means of a linear rescaling of the masses:


The new dynamical variables indeed obey


As a matter of fact, the model has a deeper dynamical symmetry. Setting


where is a global proper time, so that , i.e.,


the dynamical equations (1) can be recast as


Remarkably, the dynamics so defined are entirely parameter-free. Quadratic differential systems such as the above have attracted much attention in the mathematical literature, such as in discussions of Hilbert’s 16th problem (see e.g. Jenks (1969); Li (2003)).

Equations (6) can be formally integrated as


The amplitudes are therefore decreasing functions of . For each node , either of two things might happen:

  • Node survives asymptotically. This occurs when the integral in (7) converges in the limit. The amplitude reaches a non-zero limit , so that the mass grows exponentially as

  • Node does not survive asymptotically. This occurs when the integral in (7) diverges in the limit. This divergence is generically linear, so that the amplitude falls off to zero exponentially fast in , while the mass falls off as a double exponential in time .

The dynamics therefore drive the system to a non-trivial attractor, i.e., an extended pattern of survivors and non-survivors. This attractor depends on the whole initial mass profile (although it is independent of the overall mass scale). The formula (7) generically implies the following local constraints:

  1. each survivor is isolated (all its neighbors are non-survivors),

  2. each non-survivor has at least one survivor among its neighbors.

Conversely, every pattern obeying the above constraints is realized as an attractor of the dynamics, for some domain of initial data. This situation is therefore similar to that met in a variety of statistical-mechanical models ranging from glasses to systems with kinetic constraints. Attractors play the role of metastable states, which have been given various names, such as valleys, pure states, quasi-states or inherent structures Thouless et al. (1977); Kirkpatrick and Sherrington (1978); Stillinger and Weber (1982); Kirkpatrick and Wolynes (1987); Franz and Virasoro (2000). In all these situations the number of metastable states grows exponentially with system size as


where is the configurational entropy or complexity. This quantity is not known exactly in general, except in the one-dimensional case where it can be determined by means of a transfer-matrix approach (see Appendix A). It is relevant to mention the Edwards ensemble here, which is constructed by assigning a thermodynamical significance to the configurational entropy Edwards (1994). According to the Edwards hypothesis, all the attractors of a given ensemble (e.g. at fixed survivor density) are equally probable. This hypothesis holds generically for mean-field models, while it is weakly violated for finite-dimensional systems Franz and Virasoro (2000); Dean and Lefèvre (2001); Berg and Mehta (2001); Prados and Brey (2001); de Smedt et al. (2002); Godrèche and Luck (2005).

Iii Exact results for small systems

In this section we consider the model on small one-dimensional graphs, i.e., closed rings and open chains of  nodes. In one dimension, (1) and (6) read


for , with appropriate boundary conditions: periodic (, ) for rings and Dirichlet () for chains.

Equations (11) are very reminiscent of those defining the integrable Volterra chain. The coupling term involves the sum in the present model, whereas it involves the difference in the Volterra system. The appearance of a difference is however essential for integrability Mikhailov et al. (1987); Yamilov (2006). The present model is therefore not integrable, even in one dimension.

Our goal is to characterize the attractor reached by the dynamics on small systems of increasing sizes, as a function of the initial mass profile. This task soon becomes intractable, except on very small systems, due to the intrinsic complexity of the model. The numbers and of these attractors on rings and chains of  nodes are given in Table 2 of Appendix A. These numbers grow exponentially fast with , according to (9), with given by (86).

Ring with

The system consists of two nodes connected by two bonds. Both attractors consist of a single survivor. The dynamical equations (11) read


The difference is a conserved quantity. If , the attractor is (meaning that only node 1 survives) and its final amplitude is


and vice versa. The survivor is always the node with the larger initial mass.

In the borderline case of equal initial masses, the integrals in (7) are marginally (logarithmically) divergent. Both masses saturate to the universal limit , irrespective of their initial value.

Chain with

The two nodes are now connected by a single bond. The dynamical equations are identical to (12), up to an overall factor of two. Here too, the more massive node is the survivor.

Ring with

The attractors again consist of a single survivor. Although there is no obviously conserved quantity, the attractor can be predicted by noticing that


The sign of any difference is therefore conserved by the dynamics. In other words, the order of the masses is conserved. In particular the survivor is the node with the largest initial mass.

Chain with

The central node 2 plays a special role, so that the two attractors are and . There are two conserved quantities, and . If , the attractor is and the final amplitudes read


If , the attractor is and .

Ring with

The two attractors are the ‘diameters’ and . The alternating sum is a conserved quantity. If , the attractor is and . If , the attractor is and . The attractor is therefore always the diameter with the larger total initial mass. The individual asymptotic amplitudes cannot however be determined in general.

Chain with

The three attractors are , and . The alternating sum is a conserved quantity. This is the first case where the attractor cannot be predicted analytically in general.

Figure 1 shows the attractors reached as a function of and , for fixed . It is clear that can only be reached for , i.e., above the diagonal, whereas can only be reached for , i.e., below the diagonal. The intermediate pattern is observed in a central region near the diagonal. The form of this region can be predicted to some extent. On the horizontal axis, the transition from to takes place for . Similarly, on the vertical axis, the transition from to takes place for . The central region where is the attractor shrinks rapidly with increasing distance from the origin. This can be explained by considering the dynamics on the diagonal, i.e., in the symmetric situation where and . These symmetries are preserved by the reduced dynamics


The reduced attractor is , the full attractor is , and vanishes identically. The reduced dynamics have another conserved quantity, . The asymptotic amplitude becomes exponentially small as increases. The width of the central green region is expected to follow the same scaling law, i.e., to become exponentially narrow with distance from the origin, in agreement with our observation.

Figure 1: (Color online) Attractor on the chain with in the plane, for fixed .

Ring with

There are five attractors consisting of two survivors, obtained from each other by rotation: , , , and . There is no obviously conserved quantity, and the attractor cannot be predicted analytically in general.

Figure 2 shows the attractors as a function of and , for fixed , and . The five attractors meet at point P . If launched at P, the system is driven to the unique symmetric solution where all masses converge to the universal limit . A linear stability analysis around the latter solution reveals that its stable manifold is three-dimensional, in agreement with the observation that its intersection with the plane of the figure is the single point P.

Figure 2: (Color online) Attractor on the ring with in the plane, for fixed , and . The five attractors meet at point P.

This is the first case which manifests one of the most interesting features of the model, that of ‘winning against the odds’ Thyagu and Mehta (2011); Mehta (2013); Aich and Mehta (2014). Numerical simulations with two structureless distributions of initial masses, uniform (uni) and exponential (exp), yield the following observations:

  • The probability that the node with largest initial mass is a survivor is 0.849 (uni) or 0.937 (exp);

  • The probability that the attractor corresponds to the largest initial mass sum among the possible attractors is 0.791 (uni) or 0.891 (exp);

  • The probability that the node with smallest initial mass is a survivor is 0.018 (uni) or 0.019 (exp).

This small but non-zero probability for the smallest initial mass to survive is the beginning of the complexity associated with the phenomenon of winning against the odds. It happens essentially because more distant nodes can destroy massive intermediaries between themselves and the small masses concerned, thereby letting the latter survive ‘against the odds’. For larger sizes, the problem soon becomes intractable. First, the number of attractors grows exponentially fast with ; second, larger system sizes make for increased interaction ranges for a given node. This makes it increasingly probable to have winners against the odds, making it more and more difficult to predict attractors based only on initial mass distributions.

The above algorithmic complexity goes hand in hand with the violation of the Edwards hypothesis as well as other specific out-of-equilibrium features of the attractors, including a superexponential spatial decay of various correlations Luck and Mehta (2005). This phenomenon, put forward in zero-temperature dynamics of spin chains de Smedt et al. (2002), is reminiscent of the behaviour of a larger class of fully irreversible models, exemplified by random sequential adsorption (RSA) Evans (1993).

Iv Approximate analytical treatment

Despite the complexity referred to above, we show here that some ‘one-body observables’ can be predicted by an approximate analytical approach. Our techniques are based on the inhomogeneous mean-field theory and rely on the assumption that the statistical properties of a node only depend on its degree  Pastor-Satorras and Vespignani (2001). Such ideas have been successfully applied to a wide class of problems on complex networks (see Barrat et al. (2008); Pastor-Satorras et al. (2015) for reviews). The thermodynamic limit is implicitly taken; also, the embedding graph is replaced by an uncorrelated random graph whose nodes have probabilities  to be connected to neighbors, i.e, to have degree . In this section, we use this framework to evaluate the survival probability of a node, given its degree and/or initial mass.

iv.1 Survival probability of a typical node

We consider first the simplest observable – the survival probability of a typical node, irrespective of initial mass or degree.

From the reduced dynamical equations (6), the initial decay rate of the amplitude of node is seen to be:


If node has degree , the above expression is the sum of the initial masses of the neighboring nodes.

From a modelling point of view, this suggests a decimation process in continuous time, where nodes are removed at a rate given by their degree  at time . The initial graph is entirely defined by the probabilities for a node to have degree . Its subsequent evolution during the decimation process is characterized by its time-dependent degree distribution, i.e., by the fractions  of initial nodes which still survive at time  and have degree . The latter quantities start from


at the beginning of the process () and converge to


at the end of the process (). Indeed, as in the original model, survivors are isolated, so that their final degree is zero. The amplitude is the quantity of interest, as it represents the survival probability of a typical node. Our goal is to determine it as a function of the probabilities .

The obey the dynamical equations


The first line corresponds to the removal of a node of degree  at constant rate , while the second line describes the dynamics of its neighbors. The removal of one neighboring node adds to the fraction  at a rate proportional to while it depletes it at a rate proportional to .

The time-dependent quantity is the rate at which a random neighbor of a given node is removed at time , which is, consistent with the above, given by the average degree of a random neighbor of the node at time . This rate can be evaluated as follows (see e.g. Pastor-Satorras et al. (2015)). The probability that a node of degree has a neighbor of degree at time is, for an uncorrelated network, independent of , and given by . (The shift from to is related to the ‘inspection paradox’ in probability theory (see e.g. Feller (1966)).) The average degree of a random neighbor of the node at time  is then given by , i.e.,


Introducing the generating series


we see that obeys the partial differential equation


with initial condition


Hence it is invariant along the characteristic curves defined by


This differential equation can be integrated as




We have therefore


For infinitely long times, irrespective of , the parameter  which labels the characteristic curves converges to the limit


which we call the dynamical fugacity of the model.

We are left with the following simple expression for the survival probability of a typical node (see (19))


Furthermore, the fugacity can be shown to be implicitly given by


where the accent denotes a derivative.

We list a few quantitative predictions for important graphs/networks below:

  • Erdös-Rényi (ER) graph. This historical example of a random graph Erdös and Rényi (1959, 1960) has a Poissonian degree distribution of the form


    so that and . We obtain

  • -regular graph. In the case of a -regular graph, where all nodes have the same degree , we have . We obtain


    For the above results become

  • Geometric graph. In the case of a geometric degree distribution with parameter , i.e.,


    we have and . We obtain

  • Barabási-Albert (BA) network. The BA network, grown with a linear law of preferential attachment, has a degree distribution Barabási and Albert (1999); Barabási et al. (1999)


    with a power-law tail with exponent . We have and


    Solving (31) numerically, we obtain

  • Generalized preferential attachment (GPA) network. This is a generalization of the BA network, where the attachment probability to an existing node with degree is proportional to  Dorogovtsev et al. (2000); Krapivsky et al. (2000); Godrèche et al. (2009), with the offset representing the initial attractiveness of a node. The degree distribution


    has a power-law tail with a continuously varying exponent . We have, as expected for a tree, , irrespective of , and


    where is the Gauss hypergeometric function.

    The GPA network can be simulated efficiently by means of a redirection algorithm Krapivsky and Redner (2001, 2002). Every new node is attached either to a uniformly chosen earlier node with probability , or to the ancestor of the latter node with the redirection probability


    For (i.e., ), we have a uniform attachment rule, yielding and , so that


    For (i.e., ), the model becomes singular. The converge to , while we still have , formally. In this singular limit we have


    Finally, the BA network is recovered for (i.e., ).

Figure 3 provides a summary of the above results. The survival probability of a typical node is plotted against the mean degree , for all the above. Turning to the specific example of the GPA networks, we note that (48) and (49) provide lower bounds for the dynamical fugacity  and the survival probability respectively, which both increase to their upper bounds (50) as functions of the redirection probability . This trend is reflected in the blue arrowheads in Figure 3, where the blue square corresponds to the value for the BA network.

Figure 3: (Color online) Survival probability of a typical node, as predicted by our approximate analysis, against mean degree  for several examples of graphs and networks. For the ER and geometric graphs, the predictions for continuously varying are shown as the black (lower) and red (upper) curves. Green filled circles correspond to -regular graphs for integer . The blue square shows the prediction (44) for the BA network, while the vertical bar with arrowheads shows the range of values of for GPA networks, given by the bounds (49) and (50).

In order to test the above for a few simple cases, we have measured the survival probability of a typical node in three different geometries: the 1D chain, the 2D square lattice and the BA network. We solved the dynamical equations (6) numerically, with initial masses drawn either from a uniform (uni) or an exponential (exp) distribution. The measured values of are listed in Table 1, together with our approximate predictions (36), (38) and (44).

Geometry 1D 2D BA
0.4360 0.3755 0.6767
0.3679 0.2500 0.5530
Table 1: Survival probability of a typical node of the 1D chain, the 2D square lattice and the BA network. Comparison of numerical results for uniform (uni) and exponential (exp) mass distributions and approximate analytical results (pred).

Our numerical results suggest that the survival probability depends only weakly on the mass distribution, as long as the latter is rather structureless. While our predictions are systematically lower than the numerical observations, they do indeed reproduce the global trends rather well.

It is worth recalling here that our analysis relies on a degree-based mean-field approach. Such approximate techniques are expected to give poor results in one dimension, and to perform much better on more disordered and/or more highly connected structures. A systematic investigation of the accuracy of the present approach, including the comparison of various lattices with identical coordination numbers but different symmetries, and of various random trees or networks with the same degree distribution but different geometrical correlations, would certainly be of great interest. While such extensive numerical investigations are beyond the scope of the present, largely analytical work, we hope they will be taken up in future.

iv.2 Degree-resolved survival probability

The above degree-based mean-field analysis can be extended to quantities with a richer structure. In this section we consider the degree-resolved survival probability of a node whose initial degree  is given. The degree distribution of this special node obeys the same dynamical equations (20) as those of a typical node:


with the specific initial condition . The above dynamical equations along the same lines as in Section IV.1. The generating series


obeys the partial differential equation (23), with initial condition


Using again the method of characteristics, we get, instead of (30), the following simple behavior for the degree-resolved survival probability:


where the dynamical fugacity is given by (31). The form of this suggests a simple physical interpretation: the fugacity measures the tendency of a given node to ‘escape’ annihilation. More quantitatively, is the price per initial neighbor which a node has to pay in order to survive forever. By averaging the expression (54) over the initial degree distribution , we recover the result (30) for the survival probability of a typical node.

We now introduce the survival scale


so that our key result (54) reads


This representation makes it clear that the survival scale  corresponds to the degree of the most connected survivors.

We use this to probe the survival statistics of highly connected graphs. Here, the survival scale is large as a consequence of a large mean degree ; correspondingly, there is a decay in the survival probability of a typical node, as shown in the two examples below:

  • ER graph. Here, , and so (33) and (34) yield


    The survival scale grows almost linearly with , while the survival probability falls off as .

  • Geometric graph. Equations (40) and (41) yield


    Both the growth of the survival scale and the decay of the survival probability are slower than in the ER case. The survival scale diverges sublinearly with , with an exponent , while the survival probability decays with an exponent .

Figure 3 shows that these differing trends for the ER and geometric graphs are already evident even for low connectivity . This is because of the interesting circumstance that the behavior of the degree distribution  for relatively small degrees () determines both the growth of the survival scale  as well as the decay of the survival probability in highly connected graphs.

Assuming this regime is described by a scaling form


governed by an exponent , we obtain after some algebra




The exponents and which enter the power laws (60), (61) are always smaller than one (the value corresponding to ), when we find and . This is e.g. the case for the ER graph (see (57)).

In order to test our key prediction (54), the dependence of the survival probability of a node of the BA network on its degree was computed numerically for an exponential distribution of initial masses. Figure 4 shows a logarithmic plot of against degree . We find excellent qualitative agreement with our predictions of exponential dependence, although the measured slope, corresponding to is larger than our analytical prediction of (see (44)).

Figure 4: (Color online) Logarithmic plot of the measured degree-resolved survival probability against degree for the BA network. Straight line: least-squares fit of all the data points with slope .

iv.3 Mass-resolved survival probability

Here, we extend our analysis to the mass-resolved survival probability of a node whose initial mass is given. Since this will only enter through the reduced mass , defined as the dimensionless ratio


we denote the mass-resolved survival probability by .

Along the lines of Section IV.1, we derive the following dynamical equations for the degree distribution of the special node:


with initial condition (see (18)). To compute the rate at which a random neighbor of the special node is removed, we argue as follows. This rate is a product of the probability of finding a random neighbor of degree , given by , and the rate of removal of this node. The latter is nothing but its effective degree: while the initial masses of its typical neighbors can be taken to be , the special node has an initial mass of , so that the effective degree of this node is . The average degree of a random neighbor of the special node at time is then clearly given by , i.e.,


After some algebra, the mass-resolved survival probability reduces to:


where the mass-dependent fugacity is


This formula can be recast as


where is the fugacity given by (31). We have, consistently, for .

The mass-resolved survival probability is an increasing function of the reduced mass . We compute this explicitly for the classes of random graphs and networks considered in Section IV.1.

  • ER graph. In this case, we have



  • -regular graph. In this case, we have




    in the generic case (), whereas


    for .

  • Geometric graph. In this case, we have



  • BA and GPA networks. In these cases, numerical values of and of can be extracted from (31) and (69), using the expressions (43), (46) of the generating series .

The dependence of the survival probability of a node on its initial mass was computed numerically for the 1D chain, the 2D square lattice and the BA network, with an exponential mass distribution. Figure 5 shows plots of the measured values of against . The dashed curves show the prediction (67), rescaled so as to agree with the numerics for . In the three geometries considered, the analytical prediction reproduces the overall mass dependence of reasonably well. The observed dependence is slightly more pronounced than predicted on the 1D and 2D lattices, while the opposite holds for the BA network.

Figure 5: (Color online) Mass-resolved survival probability against reduced mass . Full curves: numerical results. Corresponding dashed curves: analytical predictions. Top to bottom near : BA network, 1D chain, 2D square lattice.

Last but by no means least, there is a striking manifestation of (super-)universality in the regime of large reduced masses, when (68) simplifies to


The mass-resolved survival probability then goes to unity according to the simple universal law


This key result in the large-mass limit is one of the strongest results in this paper: all details of the structure and embeddings of networks disappear from the survival probability of a node, leaving only a simple dependence on its mass and the mean connectivity .

iv.4 Degree and mass-resolved survival probability

Finally, our analysis can be generalized to the full degree and mass-resolved survival probability of a node whose initial degree and reduced mass are given.

The degree distribution of this special node obeys


where the rate is given by (66), and with the specific initial condition . After some algebra along the lines of the previous sections, we find:


This last result encompasses all the previous ones, including the expression (54) for the degree-resolved survival probability and the expression (67) for the mass-resolved survival probability .

The detailed numerical evaluation of degree and mass-resolved data on networks is deferred to future work, although we expect the levels of agreement to to be similar to those obtained in Sections IV.2 and IV.3. Our point of emphasis here is the simplicity of (80), which shows that the dynamical fugacity is absolutely the right variable to highlight an intrinsic universality in this problem. In showing that the result (30) is resolvable into components of degree and mass, it also points the way towards a deeper understanding of the ‘independence’ of these parameters in the survival of a node.

Finally, and no less importantly, the result (80) again manifests (super-)universality in the regime of large reduced masses (). As a consequence of (77), the degree and mass-resolved survival probability goes to unity according to the simple law


The beauty of this result (as well as its analogue (78)) lies in the fact that it is exactly what one might expect intuitively; it suggests that the probability () that a node of mass and degree might not survive, is directly proportional to its degree and inversely proportional to its mass. In everyday terms, the lighter the node, and the more well-connected it is, the more it is likely to disappear. The emergence of such startling, intuitive simplicity in an extremely complex system is a testament to an underlying elegance in this model.

V Discussion

The problem of finding even an approximate analytical solution to a model which contains multiple interactions is very challenging. In the context of predator-prey models, the (mean-field) Lotka-Volterra dynamical system and the full Volterra chain are among the rare examples which are integrable. Most other nonlinear dynamical models with competing interactions are not, and are quite simply intractable analytically.

The model inspired by black holes Luck and Mehta (2005), on which this paper is based, shows how competition between local and global interactions can give rise to non-trivial survivor patterns, and to the phenomenon referred to as ‘winning against the odds’. That is, a given mass can win out against more massive competitors in its immediate neighborhood provided that they in turn are ‘consumed’ by ever-more distant neighbors. When it was found numerically Thyagu and Mehta (2011); Mehta (2013); Aich and Mehta (2014) that such survivor distributions seemed to exhibit somewhat surprising features of universality, it was natural to ask the question: could one find the reasons for such behaviour, in the sense of characterizing these distributions at least approximately from an analytic point of view? An additional motivation was found in the work of the Barabási group Wang et al. (2013) on citation networks, where the authors put forward a universal scaling form for the ‘survival’ of a paper in terms of its citation history.

The black-hole model as defined in its original cosmological context Majumdar (2003); Majumdar et al. (2005) had an explicit time dependence due to the presence of a ubiquitous ‘fluid’, as well as a threshold below which even isolated particles did not survive. Neither of these attributes was necessary for the behaviour of most interest to us, namely the multiplicity of attractors (which in our case involve survivor distributions), and their non-trivial dependence on the initial mass profile, as a result of multiple interactions. One of our major achievements in this paper has been the construction of a much simpler model (without the unnecessary complications referred to above), which still retains its most interesting features from the point of view of statistical physics.

Once derived and established, this simple model was the basis of our investigations of universality in survivor distributions. The exact characterisation of attractors as a function of the initial data becoming rapidly impossible, we were led to think of approximate analytical techniques. Our choice of the inhomogeneous (or degree-based) mean-field theory was motivated by our emphasis on random graphs and networks in earlier numerical work Thyagu and Mehta (2011); Mehta (2013); Aich and Mehta (2014). This approach was embodied in an effective decimation process. Some of the analytical results so obtained were robustly universal, including the exponential fall-off (54) of the survival probability of a node with its degree, or the asymptotic behaviours (78) and (81) in the large-mass regime.

Our approach led us to introduce the associated concept of a dynamical fugacity, key to unlocking the reason behind the manifestation of universality in diverse survivor distributions. Physically, this signifies the tendency of a typical node in a network to escape annihilation, which we illustrate via a simple argument. Every time an agent encounters other, potentially predatory, agents it pays a price in terms of its survival probability: as the probability of each such encounter is independent of the others, the ‘cost’ to the total probability is multiplicative in terms of the number of predators encountered. The dynamical fugacity is then nothing but the cost function per encounter (i.e., per neighbor), so that the survival probability of the original agent depends exponentially on the number of its competitors. In a statistical sense, this depends only on the degree distribution of the chosen network, leading to the emergence of a universal survival probability for a given class of networks. The asymptotic survival probability of very heavy nodes becomes ‘super-universal’, in the sense of losing all dependence on different geometrical embeddings. Its form also has an appealing simplicity, as the complement of the ratio (degree/mass) of a given node; the heavier the node and the more isolated it is, the longer it is likely to survive.

In conclusion, we have used inhomogeneous mean-field theory to formulate and solve analytically an intractable problem with multiple interactions. While our analytical solutions are clearly not exact (due to the technical limitations of mean-field theory) they are nevertheless the only way to date of understanding the behaviour of the exact system. In particular, and importantly, our present analysis strongly reinforces the universality that has indeed been observed in earlier numerical simulations of this problem Thyagu and Mehta (2011); Mehta (2013); Aich and Mehta (2014). That such universal features emerge in a highly complex many-body problem with competing predatory interactions, is nothing short of remarkable.

It is a pleasure to thank Olivier Babelon and Alfred Ramani for interesting discussions.

Appendix A Numbers of attractors and complexity in one dimension

In this appendix we investigate the attractor statistics of the one-dimensional problem by means of the transfer-matrix formalism. We describe an attractor as a sequence of binary variables or spins:


From a static viewpoint, attractors are defined as patterns obeying the constraints listed in Section II. They can therefore be identified with sequences which avoid the patterns 11 and 000. The last two symbols of such a sequence may therefore be 00, 01 or 10. The numbers , and of attractors of length of each kind obey the recursion


where the transfer matrix reads


Its characteristic polynomial is , and so the Cayley-Hamilton theorem implies the recursion


Hence all the numbers grow exponentially with , in agreement with (9). The complexity reads


with being the largest eigenvalue of , i.e., the largest root of .

The total numbers of attractors on chains and on rings of nodes obey recursions derived from (85), i.e.,


with two different sets of initial conditions. The sequences and are listed in the OEIS OEIS (), respectively as entries A001608 and A000931, together with many combinatorial interpretations and references. The first few terms are listed in Table 2.

0 2 3 2 5 5 7 10 12 17 22 29
1 2 2 3 4 5 7 9 12 16 21 28
Table 2: Numbers of attractors on rings and on chains of nodes, for up to 12.
Figure 6: (Color online) Density-dependent static complexity against survivor density in the allowed range ().

The transfer-matrix approach can be generalized in order to determine the density-dependent static complexity , characterizing the exponential growth of the number of attractors with a fixed density of survivors. Introducing a static fugacity conjugate to the number of survivors, the transfer matrix becomes


whose characteristic polynomial is . The reader is referred to de Smedt et al. (2002); Godrèche and Luck (2005) for details. We obtain after some algebra


Figure 6 shows a plot of this quantity against survivor density in the allowed range (). The expression (89) reaches a maximum equal to (see (86)) when equals the mean static density of survivors:



  • Murray (1989) J. D. Murray, Mathematical Biology (Springer, Berlin, 1989).
  • Berryman (1992) A. A. Berryman, Ecology 73, 1530 (1992).
  • Szabó and Fáth (2007) G. Szabó and G. Fáth, Phys. Rep. 446, 97 (2007).
  • Majumdar (2003) A. S. Majumdar, Phys. Rev. Lett. 90, 031303 (2003).
  • Majumdar et al. (2005) A. S. Majumdar, A. Mehta, and J. M. Luck, Phys. Lett. B 607, 219 (2005).
  • Luck and Mehta (2005) J. M. Luck and A. Mehta, Eur. Phys. J. B 44, 79 (2005).
  • Thyagu and Mehta (2011) N. N. Thyagu and A. Mehta, Physica A 390, 1458 (2011).
  • Aich and Mehta (2014) S. Aich and A. Mehta, Eur. Phys. J. Special Topics 223, 2745 (2014).
  • Mehta (2013) A. Mehta, in Econophysics of Systemic Risk and Network Dynamics, edited by F. Abergel, B. K. Chakrabarti, A. Chakraborti, and A. Ghosh (Springer, Milan, 2013), New Economic Windows, pp. 141–156.
  • Mehta et al. (2005) A. Mehta, A. S. Majumdar, and J. M. Luck, in Econophysics of Wealth Distributions, edited by A. Chatterjee, B. K. Chakrabarti, and S. Yarlagadda (Springer, Milan, 2005), pp. 199–204.
  • Jenks (1969) R. D. Jenks, J. Diff. Equations 5, 497 (1969).
  • Li (2003) J. B. Li, Int. J. Bif. Chaos 13, 47 (2003).
  • Thouless et al. (1977) D. J. Thouless, P. W. Anderson, and R. G. Palmer, Phil. Mag. 35, 593 (1977).
  • Kirkpatrick and Sherrington (1978) S. Kirkpatrick and D. Sherrington, Phys. Rev. B 17, 4384 (1978).
  • Stillinger and Weber (1982) F. H. Stillinger and T. A. Weber, Phys. Rev. A 25, 978 (1982).
  • Kirkpatrick and Wolynes (1987) T. R. Kirkpatrick and P. G. Wolynes, Phys. Rev. A 35, 3072 (1987).
  • Franz and Virasoro (2000) S. Franz and M. A. Virasoro, J. Phys. A 33, 891 (2000).
  • Edwards (1994) S. F. Edwards, in Granular Matter: An Interdisciplinary Approach, edited by A. Mehta (Springer, New York, 1994).
  • Dean and Lefèvre (2001) D. S. Dean and A. Lefèvre, Phys. Rev. Lett. 86, 5639 (2001).
  • Berg and Mehta (2001) J. Berg and A. Mehta, Europhys. Lett. 56, 784 (2001).
  • Prados and Brey (2001) A. Prados and J. J. Brey, Phys. Rev. E 66, 041308 (2001).
  • de Smedt et al. (2002) G. de Smedt, C. Godrèche, and J. M. Luck, Eur. Phys. J. B 27, 363 (2002).
  • Godrèche and Luck (2005) C. Godrèche and J. M. Luck, J. Phys. Cond. Matter 17, S2573 (2005).
  • Mikhailov et al. (1987) A. V. Mikhailov, A. B. Shabat, and R. I. Yamilov, Russian Math. Surveys 42, 1 (1987).
  • Yamilov (2006) R. Yamilov, J. Phys. A 39, R541 (2006).
  • Evans (1993) J. W. Evans, Rev. Mod. Phys. 65, 1281 (1993).
  • Pastor-Satorras and Vespignani (2001) R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett. 86, 3200 (2001).
  • Barrat et al. (2008) A. Barrat, M. Barthélemy, and A. Vespignani, Dynamical Processes on Complex Networks (Cambridge University Press, Cambridge, 2008).
  • Pastor-Satorras et al. (2015) R. Pastor-Satorras, C. Castellano, P. van Mieghem, and A. Vespignani, Rev. Mod. Phys. 87, 925 (2015).
  • Feller (1966) W. Feller, An Introduction to Probability Theory and its Applications (Wiley, New York, 1966).
  • Erdös and Rényi (1959) P. Erdös and A. Rényi, Publ. Mathematicae 6, 290 (1959).
  • Erdös and Rényi (1960) P. Erdös and A. Rényi, Publ. Math. Inst. Hungar. Acad. Sci. 5, 17 (1960).
  • Barabási and Albert (1999) A. L. Barabási and R. Albert, Science 286, 509 (1999).
  • Barabási et al. (1999) A. L. Barabási, R. Albert, and H. Jeong, Physica A 272, 173 (1999).
  • Dorogovtsev et al. (2000) S. N. Dorogovtsev, J. Mendes, and A. Samukhin, Phys. Rev. Lett. 85, 4633 (2000).
  • Krapivsky et al. (2000) P. L. Krapivsky, S. Redner, and F. Leyvraz, Phys. Rev. Lett. 85, 4629 (2000).
  • Godrèche et al. (2009) C. Godrèche, H. Grandclaude, and J. M. Luck, J. Stat. Phys. 137, 1117 (2009).
  • Krapivsky and Redner (2001) P. L. Krapivsky and S. Redner, Phys. Rev. E 63, 066123 (2001).
  • Krapivsky and Redner (2002) P. L. Krapivsky and S. Redner, J. Phys. A 35, 9517 (2002).
  • Wang et al. (2013) D. Wang, C. Song, and A. L. Barabási, Science 342, 127 (2013).
  • (41) OEIS, The On-Line Encyclopedia of Integer Sequences, URL http://oeis.org.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description