Self-organization and time-stability of social hierarchies

The formation and stability of social hierarchies is a question of general relevance. Here, we propose a simple generalized theoretical model for establishing social hierarchy via pair-wise interactions between individuals and investigate its stability. In each interaction or fight, the probability of “winning” depends solely on the relative societal status of the participants, and the winner has a gain of status whereas there is an equal loss to the loser. The interactions are characterized by two parameters. The first parameter represents how much can be lost, and the second parameter represents the degree to which even a small difference of status can guarantee a win for the higher-status individual. Depending on the parameters, the resulting status distributions reach either a continuous unimodal form or lead to a totalitarian end state with one high-status individual and all other individuals having status approaching zero. However, we find that in the latter case long-lived intermediary distributions often exist, which can give the illusion of a stable society. As we show, our model allows us to make predictions consistent with animal interaction data and their evolution over a number of years. Moreover, by implementing a simple, but realistic rule that restricts interactions to sufficiently similar-status individuals, the stable or long-lived distributions acquire high-status structure corresponding to a distinct high-status class. Using household income as a proxy for societal status in human societies, we find agreement over their entire range from the low-to-middle-status parts to the characteristic high-status “tail”. We discuss how the model provides a conceptual framework for understanding the origin of social hierarchy and the factors which lead to the preservation or deterioration of the societal structure.

B Time evolution of distributions B1 Beginning from an egalitarian initial condition Fig B1 shows the time evolution of systems prepared in an "egalitarian" initial condition, in which all individuals have the same initial status of S = 1. Each row in the figure corresponds to a different value of δ, and each column to an instant in time, with time increasing from left to right across the subplots. For all values of δ, the system undergoes an early-time evolution away from the initial condition. When α = 0 (blue curves), the system arrives at a steady-state in which the distribution of status remains constant in time. The steady-state distribution is obtained within a time equal to a few multiples ofτ 1 , as expected from Eq 5 of the main text. When α > 0 (green and red curves), the distribution continues to evolve after the transient period, and approaches the totalitarian end-state at a rate that depends on δ and α. , with α = 0 (blue), α = 0.2 (green), and α = 0.5 (red), starting from an initial condition in which all individuals have status S = 1 ("egalitarian initial condition"). x-axes in the sub-plots correspond to status, S, and y-axes to the probability density, p(S). N = 10 3 , n r = 20.

B2 Beginning from a uniform initial condition
Figures showing the evolution of the model beginning from a uniform initial condition are included below, for the same values of δ and α shown in Fig. B1. In the uniform initial condition, the initial status of each individual is chosen randomly from a uniform distribution with average statusS = 1. The time evolution beginning from this initial condition is qualitatively the same as for the egalitarian initial condition. , and α = 0.5 (red), starting from an initial condition in which status is uniformly distributed ("uniform initial condition"). x-axes in the sub-plots correspond to status, S, and y-axes to the probability density, p(S). N = 10 3 , n r = 20.
C Behaviour of τ 1 and τ 2 as a function of system size, N In order to examine the behaviours of the characteristic times τ 1 and τ 2 as a function of system size, N , plots of M 2 (t) for different system sizes are shown in Fig C1. The parameters controlling the early-time behaviour of M 2 (c 1 and τ 1 ) remain constant as system size, N , is increased. The parameters controlling the long-time behaviour of M 2 increase linearly with N . For example, for α = 0.5, the values of τ 2 extracted from a fit of Eq 7 of the main text are 1.17 × 10 5 , 1.18 × 10 6 , and 1.21 × 10 7 for system sizes N = 10 2 , 10 3 , and 10 4 , respectively. Likewise, the parameter c 2 , which represents the "upper plateau" end-state value of M 2 increases linearly with N in the large-N limit (Eq 6 of the main text). D Proof that N 2 fights are required to reach end-state when δ = 1 and α = ∞ In order to explore how the long-time behaviour of the evolution of the status distribution depends on system size N , we consider the extreme scenario of a system in which the higher-status individual in each fight wins with probability p = 1, and in which the winner of each fight receives all of the loser's status, leaving the loser with zero status. This extreme scenario corresponds to δ = 1 and α = ∞. In the model presented in the main text, δ is strictly less than 1, meaning that no single individual can ever possess all of the societal status of the system. However, in the following we set δ = 1 in order to allow the system to obtain (after a finite time and for a finite system size) the end-state in which a single individual possesses all of the societal status and all other individuals have zero status. This simplifies the analysis, leading to the proof described below.
The goal of this exercise is to determine the time, τ end , required to reach the totalitarian end-state beginning from an egalitarian initial condition. To do so, we will work backwards from the totalitarian end-state, in which a single individual possesses all of the societal status.
Working backwards from the end-state, we can consider the sequence of events that were required to arrive at the end-state. The first required event, working backwards, is a fight between two individuals possessing non-zero status. Since δ = 1, one of these individuals loses and exits the fight with status equal to zero. The other individual wins the fight and receives all of the status of the loser, which, added to its pre-fight status, equals the total status of the system. Therefore, the end-state is preceded by a configuration of the system in which only two individuals possess non-zero status. The probability that these two particular individuals are selected to fight each other is the probability of selecting the first individual OR the second individual AND the probability of selecting the remaining individual: . (S1.1) Continuing to work backwards, in order to obtain the configuration of the system in which there are two individuals with non-zero status, a preceding configuration having three individuals with non-zero status is required. The probability that two of these three individuals are selected to fight is: . (S1.2) Following through with this logic, we find ρ N = N 2 2 N (N −1) . Now, the expected number of fights that will elapse between the required configuration k and the subsequent required configuration k − 1 in this sequence is ρ −1 k . Therefore, the expected number of fights to reach the end-state is: The sum in Eq. S1.3 approaches 1 for large N , such that τ end ≈ N 2 . The time to reach the end-state when δ = 1 and α = ∞ is therefore τ end ≈ 2τ end /N = 2N , for large N .
This proof is included as a demonstration of the configurational reasons why the long-time (approach to the end-state) evolution of the model dynamics increases in proportion to the system size N .  Three distinct stages in the evolution of M 2 (t) can be identified, as indicated on the log-log plot of M 2 (t) shown in E1a. First (stage 1), there is a rapid change away from the initial condition into a distribution similar in shape to the α = 0 (true steady-state) distribution. This distribution shape is essentially maintained during stage 2, such that M 2 (t) changes only very slowly with time. For example the curve corresponding to α = 0.2 (blue) in Fig E1a remains in stage 2 over the time scale of the simulation (to demonstrate this, distributions corresponding to specific points in time along this curve are shown in Fig E2). The duration of stage 2 decreases with increasing α, as can be seen in the curve corresponding to α = 0.4 (green) in Fig  E1a (distributions corresponding to specific points in time along this curve are shown in Fig E3). When α is increased further, stage 2 essentially disappears (red curve in Fig E1a -distributions corresponding to specific points along this curve are shown in Fig E5). Transitioning from stage 2 to stage 3, the distribution accumulates many individuals with very low statuses, and several higher-status individuals emerge. Stage 3 then consists of the slow evolution of this highly skewed distribution. A particular characteristic of stage 3 is that, for fights involving at least one individual who still retains sizeable status (e.g., S 1 >S), the probability p → 1. This is illustrated in Fig E1b, in which p is averaged over the preceding 100 fights, considering only those fights in which S 1 >S, giving the quantity p S1>S . The fact that p S1>S → 1 signifies that the higher-status individual in any given fight (excluding fights in which neither individual has greater than average status) is virtually guaranteed to win. For finite values of α > 0, in order for this to occur, there must be a large relative status between the two individuals. Stage 3 is therefore characterized by a distribution in which the ratio of S 2 /S 1 for any two randomly-selected individuals is small enough such that p S1>S → 1 for fixed α (as per Eq 1 of the main text). The large-time decrease in p S1>S (e.g. for the red curve corresponding to α = 0.75 in Fig E1b, at ln[t] = 8 to ln[t] = 11) occurs as individuals withS < S < S H begin to experience decreases in status. Here, S H is the status of the highest-status individual. This is illustrated in Fig E4, which plots the sum of statuses of all individuals withS < S < S H as a function of time. The transition from stage 3 to the totalitarian end-state involves the decrease in the sum of statuses shown in Fig E4 following the maximum at ln[t] ≈ 7. This shows that in its evolution toward the end-state, the distribution must pass through a series of configurations (stage 3) in which there are a number of individuals with higher-than-average statuses, these individuals having extracted their large statuses from the mass of very low status individuals that accumulates in the transition from stage 2 to stage 3. Only after stage 3 has been obtained are the conditions present in which fights between high status individuals can lead to the emergence, in the end-state, of a single individual with status approaching the total status of the system.  A plot showing the same quantities as in Fig E1, but where α is fixed and δ is varied, is included below for comparison. As can be seen, as δ is increased, the duration time of stage 2 decreases. Consider a large (but finite) system of N individuals. When α > 0, beginning from an egalitarian initial condition, the system evolves such that there are relatively few individuals with large status and relatively many individuals with small status (see section E). Since the average status of the system is constant, we can use it as a reference point to define "large status" to mean greater than average status, and "small status" to mean less than average status. The subsequent dynamics can be understood in terms of two timescales. The first timescale, λ 1 , sets the time (measured in number of fights) that transpires, on average, between fights involving the same two particular individuals. Since pairs of individuals are picked randomly from among the population, λ 1 is determined by the system size N and remains unchanged as the system evolves.
The second timescale, λ 2 , relates to the time required for an individual with less than average status to achieve an upset victory over an individual with greater than average status. Once the system has evolved to a point in which there are few individuals with high statuses and many individuals with low statuses, such an upset is the only way in which an individual with lower than average status can obtain greater than average status. λ 2 is determined by the average probability, p upset , with which an upset occurs (the smaller p upset , the larger λ 2 ). The probability for such an upset to occur (in any particular fight) is less than 50%, by Eq 1 of the main text. Therefore, fights that involve one individual with greater than average status and one individual with less than average status will generally result in a transfer of status from the lower-status individual to the higher-status individual. This creates a net transfer of status from individuals with lower than average status to individuals with greater than average status. One outcome of this net transfer of status is a reduction in p upset . Consequently, the timescale λ 2 must increase with time (number of fights).
Since λ 1 remains constant whereas λ 2 increases as the system evolves, eventually λ 1 dominates the dynamics such that any two particular individuals fight one another much more frequently than any one particular individual with lower than average status achieves an upset that transforms it into an individual with greater than average status. Since fights between any two particular high-status individuals occur with higher frequency than fights in which a low-status individual achieves an upset over a high-status individual, there is a net flow of the losers of fights between two high-status individuals into the ranks of those with lower than average status.
The rate at which this net flow of losers occurs increases as λ 2 increases, and after an infinite number of fights (t → ∞), only a single individual with greater than average status remains.

G Technical aspects regarding τ 2 extraction and errors
For small values of α, τ 2 diverges, such that very long simulations are required in order to evaluate τ 2 . In order to simplify matters, we consider a linear expansion of Eq 7 of the main text in the region where τ 1 t τ 2 : A linear fit to M 2 (t) data in the τ 1 t τ 2 can then be made, with intercept c 1 and slope s = (c 2 −c 1 )/τ 2 . This allows τ 2 to be calculated, since c 2 is known from Eq 6 of the main text.
An example is shown in Fig G1a. Here, a set of 10 realizations of the simulation were performed. The plot contains the 10 different M 2 (t) curves overlayed one top of each other, with a linear fit (red line) to the combined data for all 10 curves. Fig G1b shows the distribution of residuals. The residuals are skewed, such that large increases in M 2 are more likely than large decreases.
Propagation of errors shows that ∆ ln[N/τ 2 ] ≈ ∆s/s, assuming that ∆s ∆c 1 (justified below). This means that as α is decreased and s approaches 0, it eventually becomes impossible to obtain a meaningful fitted value of τ 2 . Following this reasoning, we include points on the plot of ln[N/τ 2 ] vs. α shown in Fig 5a  of the main text for values of α for which ∆s/s < 1.
Since M 2 (t) has a skewed (non-Normal) distribution of residuals and contains time-correlations, it is difficult to obtain a meaningful value of ∆s from a regression. We therefore estimate ∆s using uncorrelated, normally-distributed synthetic data with a width that is representative of the fluctuations in M 2 (t). This is done by generating random data from a Normal distribution with width equal to the average (over all time points) of the standard deviation taken, at each time t, across the 10 realizations of the simulation. An example of this synthetic data is shown in Fig G1b (black curve). A least squares fit provides a correct determination of the error on the slope and intercept of this synthetic data, which we use as estimates for ∆s and ∆c 1 . In this way, we see that ∆s ∆c 1 .
This procedure provides an estimate of the (small α) limit beyond which points cannot be added to the plot of Fig 5a of the main text (obtaining further points would require longer simulation times that are beyond the scope of this study). The procedure does not provide a reasonable estimate of error bars for

H Equi-σ distributions
Essentially identical distributions can be produced (within a particular simulation time) for different combinations of the parameters δ and α. This is shown in Fig H1, which plots status distributions generated for three different points in δ − α parameter-space. The inset to Fig H1 shows the distributions with a logarithmic scale on the y-axis, in order to allow closer inspection of the tails of the distributions. To compare the shapes of the distributions in Fig H1, in addition to visual inspection of the plotted curves, we use the time-and-realization-averaged standard deviation, σ of the status distributions. This is defined as follows: first, we average the standard deviation of the status distributions over 500 realizations of the simulation, giving us σ(t). Then, we average σ(t) over a range of time (500 ≤ t ≤ 1000) during which σ(t) is approximately constant corresponding to the long-lived state. This gives us σ . Repeating this procedure for the skewness of the status distributions gives us µ 3 /σ 3 , which we also use to compare the shapes of the distributions in Fig H1.  The value of σ for the three distributions is the same, within error (caption of Fig H1). The error value used here is the standard deviation of σ(t) taken over the range of time 500 ≤ t ≤ 1000. Likewise, µ 3 /σ 3 is the same, within error, for the three different distributions. There is, however, one subtle difference between the distributions that can be seen in Fig H1: the value of p(S) at S = 0 becomes systematically larger going from smaller (larger) to larger (smaller) α (δ). Except for this caveat, Fig H1 shows the existence of sets of (δ, α) points producing essentially identical distributions.

I Determining the location of the transition between regions II
and III in Fig 6 of the main text I1 Using a criterion based on the slope of M 2 (t) The location of this transition between regions II and III of Fig 6 of the main text depends on the time, τ obs , over which the system is observed (simulated). A simple criterion equates the onset of runaway with the appearance of a positive slope in the long-time portion of M 2 (t). Whether such a slope is observed in M 2 (t) depends on the size of the fluctuations in the simulation data. We ran n r = 25 realizations of the simulation for system size N = 1000. In order to determine if the system is in the long-lived state for a particular value of α, we first fit a straight line to M 2 (t) over a range [t 0 , τ obs ], where t 0 τ 1 such that M (t 0 ) is on the M 2 (t) plateau that follows the initial transient evolution away from the initial (egalitarian) distribution. The transition is considered to have occurred at the value of α for which the long-time part of the M 2 (t) curve increases by 2.5% over the observation (simulation) time τ obs . For the choice of N and n r used, this amount of increase corresponds to approximately 5 times the standard deviation of the residual R(t) =M 2 (t)−M 2 (t), whereM 2 (t) is the linear fit, and M 2 (t) is averaged over the n r = 25 realizations. Therefore, the distribution is considered to be long-lived over the interval [t 0 , τ obs ] as long as m(τ obs − t 0 ) < 0.025b, where m is the slope and b the intercept of the linear fit.

I2 Using the Arrhenius relation between τ 2 and α
The criterion (see above, section I1) of a 2.5% increase in M 2 (t), used to determine the location of the transition between the long-lived state and runaway, can be written in terms of Eq 7 of the main text as: where the approximation comes from the assumption that τ obs τ 1 . Rearranging Eq. S1.5 in terms of τ obs gives the following relationship: (S1.6) Now, we fix τ obs and adjust α until the relationship in Eq. S1.6 is satisfied. Inserting the relationships α b = 0.53δ −1.21 (Fig 5b of the main text) and N f 0 = δ 1.3 (lower inset of Fig 5b of the main text) into Eq 8 of the main text and rearranging, we find: where α t is the location of the transition (as defined by the criterion stated above) and τ 2 (α t ) is the value of τ 2 when α = α t .
A factor of c 1 appears twice in Eq. S1.7. While c 1 is in fact a function of δ and α, it changes slowly with α (Fig 4d of the main text). Therefore, we set c 1 to be equal to its value when α = 0, such that c 1 = δ/(1 − δ) (see Eq 3 of the main text). With these assumptions in place, Eq. S1.7 produces the solid black (τ obs = 10 4 ) and red (τ obs = 10 3 ) curves in Fig 6 of The N -dependence of α t is contained within the denominator of the argument of the logarithm in the denominator of Eq. S1.7. Fig I1 shows a plot of this part of Eq. S1.7 as a function of N . The plot becomes effectively constant for large N , showing that the N -dependence of α t vanishes for increasing N such that the location of the transition between the long-lived state and runaway effectively does not depend on system size.

I4 Alternative location of the transition between regions II and III according to Arrhenius relation
In some physical systems governed by an Arrhenius relation, the location of the operational transition between two regimes is taken by the researcher to be at the value of the control parameter (usually temperature, in physical systems) for which the transition rate τ is equal to the measurement time τ m . This is the case in the blocking transition of superparamagnetism, for example [1].
We could use this reasoning to locate the transition between the long-lived state and runaway in our model at a value of α = α * t for which τ 2 = τ obs . Re-arranging Eq. 8 of the main text and substituting in the relationships for α b and N f 0 , we find: (S1.8) In a physical experiment, the criterion for setting the location of the operational transition is primarily determined by the measurement type. For some measurement types, it may be appropriate or natural to set the location of the transition at a value of the control parameter analogous to α * t . In our system, this would correspond to identifying the transition between the long-lived state and runaway at a value of α for which M 2 = M 2 (t = τ 2 ) ≈ c 2 (1 − 1/e). This represents a significant evolution of the system toward the totalitarian end-state, which is not suitable for our purpose. Rather, our goal is to determine the value of α for which, on the scale of the observation time considered, the long-lived state plateau value of M 2 is lost and a noticeable increase in M 2 (t) is observed. We therefore place the location of the transition at α t .