Dynamics of ranking

Virtually anything can be and is ranked; people, institutions, countries, words, genes. Rankings reduce complex systems to ordered lists, reflecting the ability of their elements to perform relevant functions, and are being used from socioeconomic policy to knowledge extraction. A century of research has found regularities when temporal rank data is aggregated. Far less is known, however, about how rankings change in time. Here we explore the dynamics of 30 rankings in natural, social, economic, and infrastructural systems, comprising millions of elements and timescales from minutes to centuries. We find that the flux of new elements determines the stability of a ranking: for high flux only the top of the list is stable, otherwise top and bottom are equally stable. We show that two basic mechanisms — displacement and replacement of elements — capture empirical ranking dynamics. The model uncovers two regimes of behavior; fast and large rank changes, or slow diffusion. Our results indicate that the balance between robustness and adaptability in ranked systems might be governed by simple random processes irrespective of system details.


Introduction
Rankings are everywhere. From country development indices, academic indicators and candidate poll numbers to music charts and sports scoreboards, rankings are key to how humans measure and make sense of the world [1,2]. The ubiquity of rankings stems from the generality of their definition: Reducing the (often high-dimensional) complexity of a system to a few or even a single measurable quantity of interest [3,4], dubbed score, leads to an ordered list where elements are ranked, typically from highest to lowest score. Rankings are, in this sense, a proxy of relevance or fitness to perform a function in the Ranking lists typically have a fixed size N 0 (e.g., the Top 100 universities [51], the Fortune 500 companies [52]), so elements may enter or leave the list at any of the T observations t = 0, . . . , T − 1, allowing us to measure the flux of elements across rank boundaries [42,46,47] (for the observed values of N 0 , T see SI Table S1). We introduce two time-dependent measures of flux: the rank turnover o t = N t /N 0 , representing the number N t of elements ever seen in the ranking list until time t relative to the list size N 0 , and the rank flux F t , representing the probability that an element enters or leaves the ranking list at time t. Rank turnover is a monotonic increasing function indicating how fast new elements reach the list (Fig. 1b left; all datasets in SI Fig. S5). In turn, flux shows a striking stationarity in time despite differences in temporal scales and potential shocks to the system (SI Fig. S3). By averaging over time, the mean turnover rateȯ = (o T −1 − o 0 )/(T − 1) and the mean flux F = F t turn out to be highly correlated quantities that uncover a continuum of ranking lists (Fig. 1b right; values in SI Table S2). In The measures of rank turnover and flux reveal regularities in the stability of ranking processes [43,53].
We follow the time series of the rank R t occupied by a given element at time t [44] (Fig. 1c; all datasets in SI Fig. S2). In most systems, highly ranked elements like Harvard University and the English word 'the' never change position, showcasing the correspondence between rank stability and notions of relevance like academic prestige [7,27], grammatical function [17,50], and underlying network structure [53]. As we go down the ranking list of open systems, rank trajectories increasingly fluctuate in time.
In the least open systems where turnover and flux are low, however, low ranked elements are also stable. In the ranking of British cities by population, for example, both Birmingham and Nairn remain the most and least populated local authority areas throughout the 20th century [54]. These findings uncover a more fine-grained sense of rank stability: most systems have a stable top ranking, but only the least open systems feature stable bottom ranks as well. The rank change C, measured as the average probability that element at rank R changes between times t − 1 and t, varies between an approximately monotonic increasing function of R for open systems to a symmetric shape as systems become less open ( Fig. 1d; all datasets in SI Fig. S6).
Since the stability patterns of an empirical ranking list (as measured by rank change C) can be systematically connected to the amount of elements flowing into and out of the list, we build a model of rank dynamics based solely on simple generative mechanisms of flux (Fig. 2). Without assuming system-specific features of elements or their interactions, there are at least two ways to implement flux in rank space. Smooth (but arbitrarily large) changes in the score of an element might make it larger or smaller than other scores, causing elements to move across ranks (the way some scientists gather more citations than others [27], or how population size fluctuates due to historical events [44]). Regardless of score, elements might also disappear from the list and be replaced by new elements: young athletes enter competitions while old ones retire [36]; new words replace anachronisms due to cultural shifts [50].
We implement random mechanisms of displacement and replacement in a simple model by considering a synthetic ranking list of length N 0 embedded within a larger system of size N ≥ N 0 . At each time step of length ∆t = 1/N , a randomly chosen element moves to a randomly selected rank with probability We solve the model analytically by introducing the displacement probability P x,t that an element with rank r = R/N gets displaced to rank x = X/N after a time t ( Fig. 2b; uppercase/lowercase symbols denote integer/normalized ranks). Since for small ∆x = 1/N the probability that at time t an element has not yet been replaced is e −νt , we have Here, L t = (1 − e −τ t )/N is the (rank-independent) probability that up until time t an element gets selected and jumps to any other rank. The length of jumps is uniformly distributed, so they can be thought of as a Lévy random walk with step length exponent 0 [55] (full derivation in SI Section S4).
The probability D x,t = D(x, t)∆x that the element in rank r gets displaced to rank x after a time t (due to Lévy walks of other elements) follows approximately the diffusion-like equation where α = τ /N . Probability Px,t that element in rank r = R/N moves to x = X/N after time t (uppercase/lowercase symbols are integer/normalized ranks). Elements not replaced diffuse around x = r (with probability Dx,t) or perform Lévy walks [55] (with probability Lt). Eq. (2)  in allele populations [56,57]. The solution D(x, t) of Eq. (2) is well approximated by a decaying Gaussian distribution with mean r and standard deviation 2αr(1 − r)t, i.e. a diffusion kernel (Fig. 2b). Overall, local displacement makes elements slowly diffuse around their initial rank, while Lévy walks and the replacement dynamics reduces exponentially the probability that old elements remain in the ranking list.
An explicit expression for the displacement probability P x,t allows us to derive the mean flux and the mean turnover rateȯ where p = N 0 /N is the length of the ranking list relative to system size (see SI Section S4). In order to fit the model to each empirical ranking list, we obtain N 0 from the data and approximate N = N T −1 as the number of distinct elements ever seen in the list during the observation period T , thus fixing p (values for all datasets in SI Tables S1-S2). The remaining free parameters τ and ν (regulating the mechanisms of displacement and replacement) come from numerically solving Eqs. introduce a small bias in the estimation of τ (SI Fig. S18). Despite this bias, the simple generative   SI Table S1 and SI Section S5). Online social systems have the largest rates, followed by sports and languages (SI Table S2). (Bottom) Parameters τr and νr for T eff = T /k subsampled observations (all datasets in SI Fig. S22). By subsampling ranking dynamics, systems move downwards along the universal curve while keeping a constant replacement rate. (c)-(e) Average probability that an element changes rank by Lévy walk (W levy ), diffusion (W diff ), or is replaced (W repl ) between consecutive observations in the data. Probabilities shown both for selected datasets (dots), and for the model moving along the curve τrνr = 1 with the same p andȯ as the data (lines) (for rest of systems see SI Fig. S17). The simulated system in (e) is the model itself for given values of τ , ν, and p (shown in plot). The model reveals a crossover in real-world ranking lists between a regime dominated by Lévy walks (b) to one driven by diffusion (c). Although not seen in data, the model also predicts a third regime driven by replacement (d).
mechanisms of flux in the model are enough to recover the behavior of ranking lists as quantified by P x,t and C ( Fig. 2d and SI Fig. S19): When rank flux is low, both the top and bottom of the list are similarly stable and rank dynamics is mostly driven by an interplay between Lévy walks and diffusion.
As systems become more open, however, this symmetry gets broken due to a growing flux of elements at the bottom of the ranking list (see SI Fig. S4). Regardless of whether we rank people or animals, words or countries, the pattern of stability across a ranking list is accurately emulated by random mechanisms of flux that disregard the microscopic details of the individual system.
The characterization of flux in ranking lists with mechanisms of displacement and replacement of elements reveals regimes of dynamical behavior that are not apparent from the data alone (Fig. 3). By rescaling the fitted parameter values of the model as most open ranking lists (F,ȯ > 0) are predicted to follow the universal curve which suggests that ranking dynamics are regulated by a single effective parameter [ Fig. 3a; derivation in SI Section S5; for a discussion of the role of fluctuations on the validity of Eq. (6) see SI Figs. S18-S20]. Even if, potentially, displacement and replacement could appear in any relative quantity, adjusting the model to observations of rank flux and turnover (Fig. 1b) leads to an inverse relationship between parameters regulating their generative mechanisms. Real-world ranking lists lie in a spectrum where their dynamics is either mainly driven by score changes that displace elements in rank (high τ r and low ν r , like for GitHub software repositories [58] ranked by daily popularity), or by birth-death processes triggering element replacement (low τ r and high ν r , like for Fortune 500 companies [52] ranked by yearly revenue). While the symmetry (or lack thereof) in rank change C may seemingly imply two distinct classes of systems (see Fig. 1d and Fig. 2d), Eq. (6) reveals the existence of a continuum of open ranking lists, which can be captured by a single model with a single effective parameter.
Data on empirical ranking lists is constrained by the average time length between recorded observations, which varies from minutes to years depending on the source and intended use of the rankings ( for all datasets is listed in SI Table S1). We explore such scoping effect by subsampling data every k observations (for details see SI Section S5 and SI Figs. S21-S22). Longer times between snapshots of the ranking list lead to an increase in rank flux, turnover, and fitted parameters, such that the rate of element replacement ν/k stays roughly constant (Fig. 3b top). A conserved replacement probability per unit time, robust to changes in sampling rate, is yet another measure of rank stability: online social systems exchange elements frequently (e.g., the ranking lists by daily popularity of both GitHub software repositories and of online readers of the British newspaper The Guardian [59]), followed by sports, while languages are the most stable (values for k = 1 in SI Table S2).
The universal curve in Eq. (6) (Fig. 3c). Here, long-range rank changes take elements in and out of a short ranking list within a big system (low p), thus generating large mean flux F (see SI Table S2). Most datasets, like the yearly rankings of scientists by citations in American Physical Society journals [27,60] and of countries by economic complexity [34,61], belong instead to a diffusion regime with W diff > W levy , W repl (Fig. 3d). In this regime, a local, diffusive rank dynamics is the result of elements smoothly changing their scores and overcoming their neighbors in rank space. Under subsampling, ranking lists move downwards along the universal curve, going from a state with a certain number of Lévy walks to one more driven by diffusion (Fig. 3b bottom; all datasets in SI Fig. S22).
The model also predicts a third regime dominated by replacement (W repl > W levy , W diff ; Fig. 3e), where elements are more likely to disappear than change rank. Such ranking lists replace most constituents from one observation to the next, forming a highly fluctuating regime that we do not observe in empirical data. To showcase the crossover between regimes, we simulate the model along the universal curve of Eq. (6) while keeping p andȯ fixed in Eq. (5) (lines in Fig. 3c-e). These curves show how close systems are to a change of regime, i.e. from one dominated by Lévy walks to one driven by diffusion.
When a ranking list is close to a regime boundary, external shocks (amounting to variations in parameters τ and ν) may change the main mechanism behind rank dynamics, thus affecting the overall stability of the system.

Discussion
Ranking lists reduce the elements of high-dimensional complex systems into ordered values of a summary statistic, allowing us to compare seemingly disparate phenomena in nature and society [2,42]. The diversity of their components (people, animals, words, institutions, countries) stands in contrast with the statistical regularity of score-rank distributions when aggregated over time [9,37] [56,57], both alongside a relatively low rate of element replacement independent from the frequency at which the ranking list is measured.
Our results suggest that, even though score distributions differ across systems depending on what type of elements and interactions they have (SI Fig. S1), ranking lists have similar stability features.
What are the underlying properties of the system that enhance this similarity? An extension of our model explicitly considering the links between score and rank may help further understand the experimental evidence in this area, like the recent observation that the stability of crowdsourced rankings depends on the magnitude difference between quality scores [62]. It is also interesting to consider the observed deviations from the predictions of our model, even at the level of ranks. Rank flux for languages is not constant but decreases over time (SI Fig. S3), arguably due to the long observation period (over three centuries; see SI Table S1), variations in word sampling across decades, or even cognitive distortions at the societal scale [63]. The rank-dependence of flux for very open systems (SI Fig. S4) and the slow decay of inertia with long times we observe in most datasets (SI Fig. S7) might be better reproduced by a nonuniform sampling of elements in the mechanisms of displacement and replacement of the model. Finally, deviations in the data (indicating a departure from the assumptions of randomness and stationarity built into our model) could be used to detect shocks to the system larger than expected statistical fluctuations, such as the sudden increase in rank flux of Fortune 500 companies during financial crises (SI Fig. S3).
A more nuanced understanding of the generic features of ranking dynamics might help us limit resource exhaustion in competitive environments, such as information overload in online social platforms and prestige biases in scientific publishing [64], via better algorithmic rating tools [65]. The observation of a systematic interplay between "slower" and "faster" ranking dynamics [66] (see SI Section S6) can be refined by exploring the relationship between ranking lists associated to the same system, or by incorporating networked interactions that lead to macroscopic ordering [44,67], which may provide a deeper understanding of network centrality measures based on ranking [68]. Given that rankings often mediate access to resources via policy, similar mechanisms to those explored here may play a role in finding better ways to avoid social and economic disparity. In general, a better understanding of rank dynamics is promising for regulating systems by adjusting their temporal heterogeneity.

Data availability
For data availability see SI Section S2. Non-public data is available from the authors upon reasonable request.

Competing interest statement
A.-L.B. is founder of Foodome, ScipherMedicine, and Datapolis, companies that explore the role of networks in health and urban environments. G.I. is founder of Predify, a data science consulting startup in Mexico. C.P. and C.G. declare no competing interest.

Supplementary Information for
Dynamics of ranking

S1 Summary of notation
We consider a ranking list at times t = 0, . . . , T − 1, that is, an ordered set of N 0 elements (with N 0 constant in time) where each element i in the list at time t has a score s i (t), and elements are ordered across the list with decreasing score. Thus, elements and their order may change throughout time. The most important element (with the largest score) at time t has rank R t = 1; the least important element (with the lowest score) has rank R t = N 0 , and the ranking is the particular order of elements and scores across ranks R = 1, . . . , N 0 , which changes in time as elements vary their scores. We consider the case where we do not have access to the scores at all times, but only in a discrete set of T observations, separated in average by a real time interval (measured in days) 1 . We also consider the case where we may not have empirical data on the scores of all elements at all times, either because elements enter/leave the ranking at some point in time, or because score data is unavailable. Thus, there are N t ≥ N 0 distinct elements that have ever been in the ranking up to (and including) time t. If N t = N 0 for all t the ranking list is closed, since elements do not leave or enter the ranking list, and if N t > N 0 the ranking list is open. In what follows we characterize the temporal variability of ranking lists across many systems by analysing the flow of elements in and out of the ranking list as time goes by.

S2 Data description
We analyse 30 datasets in a wide range of structural and temporal scales, comprising several definitions of elements and scores. Table S1 lists the observed system size N T −1 (number of distinct elements ever seen in the ranking), ranking list size N 0 , number of observations T , and real time interval for all datasets considered. Table S1 also includes a system classification based on the nature of the elements in the system, and the corresponding definitions of elements and scores for each system. Social systems reflect human interactions at the individual and organizational levels. Language datasets show how word usage has changed across centuries. Economic rankings illustrate value at different scales. Systems in the infrastructure category are specific to public transport and city populations. Nature datasets gather information from biological and geological phenomena. Finally, rankings in sports capture the relative performance of players and teams according to sets of rules. We now describe the datasets in more detail.
In their original state some datasets are not homogeneous in time and size, since they do not have the same real time interval between any times t and t + 1 or the same number of elements per observation (i.e. they have a variable ranking list size). To consistently analyse all datasets, we crop the data to obtain roughly homogeneous time intervals and have a constant ranking list size across observations, while trying to retain as large T and N 0 as possible.

S2.1 Society
GitHub repositories. GitHub [1] is perhaps the most popular web-based version control repository, mostly used for source code. This dataset contains daily rankings of repositories, based on the number of users watching each project, from April 1, 2012 to December 30, 2014. 1 The time unit is arbitrary and may be anything other than days.

Dataset Element
Score

Sports
Chess players (male) [22] person Elo rating [23] 16568 13500 46 30.44 Chess players (female) [22] person Elo rating [23] 16539 12681 46 30.44 Poker players [24] person GPI score [25] 9799 1795 221 7.04 Tennis players [26] person ATP points [27] 4793 1600 400 7.00 Golf players [28] person OWGR points [29] 3632 1150 768 7.00 Football scorers [30] person # goals [ Table S1. Datasets used in this study. Characteristics of the available datasets, including the observed system size NT −1, ranking list size N0, number of observations T , and real time interval between observations . The table includes a system type based on the nature of the elements in the ranking list, as well as the corresponding definitions of elements and scores for each system, and references regarding the dataset and score system.
The Guardian readers. The Guardian [2] is a British national daily newspaper, publishing online articles (news and opinion pieces) in diverse subjects. Users registered to the website can post comments to some of the articles, which readers may 'recommend' (by clicking a 'like' button akin to those of social network sites). Focusing on the period from November 1, 2011 till May 1, 2012, we crawl articles appearing in three sections of The Guardian ('politics', 'sport', and 'comment is free'), and rank users daily by the number of comments they write (denoted 'comm'), as well as by the average number of recommends their comments receive (denoted 'recc').
Enron emails. The Federal Energy Regulatory Commission, during its investigation of the company Enron, made public about half a million emails from roughly 150 users, mostly senior managers. A cleaned, current version of the dataset is available online [3]. From this dataset we rank users (email accounts) by the number of emails sent, on a weekly basis.
Scientists. We use yearly citation data from journals of the American Physical Society to rank the most cited scientists by number of citations [4,5].
Universities. ShangaiRanking Consultancy is an independent organization dedicated to higher edu-

S2.2 Languages
We use a subset of the publicly available Google Books Ngram dataset [8,9], the result of the digitalization and conversion of millions of books in several languages (about 4% of all published books until 2009).
From this data we extract the frequency of words by year for Russian, Spanish, German, French, Italian, and English. The original dataset is case sensitive, so we merge words with different cases. We have already analyzed the temporal variability of ranks at several scales for these languages in previous studies [10,11].

S2.3 Economics
Companies. Since 1955, Fortune magazine has compiled a yearly dataset of the top 500 corporations in the world based on yearly revenue. Here we use the freely accessible archives of 1955-2005 [12].
Countries. MIT's Observatory of Economic Complexity [13] has compiled a yearly ranking dataset of the economic complexity of countries between 1962 and 2015 [14,15], where a complexity score reflects the diversity of exports of countries around the globe.

S2.4 Infrastructure
Cities. To obtain ranking datasets of cities by population, we use data from previous studies of Russia (RU) [16] and Great Britain (GB) [18].

S2.5 Nature
Hyenas. A 23-year field study has monitored the social relationships of a spotted hyena population in Kenya [19]. Using this dataset, we rank individual hyenas according to the sum of an association index describing the strength of the relationships with the rest of the hyenas over a particular year.
Regions JP. The Japan University Network Earthquake Catalog openly records earthquake events happening in Japan (JP) [20]. Following a previous study [21], we only consider earthquakes with Richter magnitude larger than 2.0, happening between July 1, 1985 and December 31, 1998, taking main-shocks and after-shocks as separate events. We rank administrative regions in Japan by both the monthly number of earthquakes happening in a region (denoted 'quakes'), and the average earthquake magnitude in a month (denoted 'quake mag').

S2.6 Sports
Chess players. Every month, the Fédération Internationale desÉchecs (FIDE, International Chess Federation) ranks top male and female players [22] based on the Elo rating system [23].
Poker players. Weekly rankings of Poker players are listed online [24]. Rankings are based on the Global Poker Index (GPI) score [25], which takes into account players who participated in tournaments in the last 18 months, including an aging factor and the scores obtained in different tournaments.
Tennis players. The Association of Tennis Professionals (ATP) [26] has an intricate ranking system based on ATP points [27], which is updated every week to rank male players in Singles.
Golf players. The Official World Golf Ranking (OWGR) publishes weekly lists of player rankings, accumulating points achieved in tournaments in the previous two years [28]. To calculate OWGR points, the points achieved by each player on all tournaments during these two years are averaged over the number of tournaments involved [29].
Football scorers. The Football World Rankings publishes a weekly list of top scorers [30], considering the number of goals scored by each player in the previous year [31].

NASCAR drivers.
We use an open dataset containing the ranking of drivers of the National Association for Stock Car Auto Racing (NASCAR), both for the Busch Series and the Winston Cup [32].
National football teams. The Fédération Internationale de Football Association (FIFA, International Federation of Association Football) [33] publishes a monthly list of national teams ranked according to the points obtained in the previous four years, using several criteria to determine FIFA points [34]. Data for female teams has a lower temporal scale, so we restrict our analysis to male teams.

S2.7 Acknowledgments and data availability
We acknowledge José A. Morales and Sergio Sánchez for data handling in the initial stages of the project.
We are grateful for data provision to Syed Haque (GitHub repositories and Enron emails), Raj Kumar Football scorers  between times t, we normalize the score of element i as the fitness with s m (t) the maximum score in the ranking at time t. We define P f as the probability that a randomly selected element i has fitness f i at time t, and P m as the maximum value of P f at time t. The normalized fitness distribution P f /P m for each time t = 0, . . . , T − 1 varies greatly between systems, having either a broad or narrow functional form (Fig. S1). In most ranking lists, P f is relatively constant in time.

S3.2 Rank dynamics
We consider the dynamics of rank R t = 1, . . . , N 0 of a given element as a function of time t. In most open ranking lists (N t > N 0 ), the time series R t at the top of the ranking (R t N 0 ) tend to be stable, while R t fluctuates more as we go down in ranking (middle rows in Fig. S2). If at some time t an element 0.0 0.5 The Guardian readers (recc) The Guardian readers (comm) Football scorers is not in the ranking, its rank R t is not observable from empirical data (because we do not know the corresponding score). For some open ranking lists even the top of the ranking is not stable, since the time scale T is larger than the typical lifetime of an element in the system (see, e.g., upper rows in Fig. S2, or datasets where people or companies have not been active during the whole observation period). In closed ranking lists (N t = N 0 ), both the top and the bottom of the ranking are stable, while the middle part of the ranking shows more variation in R t (lower rows in Fig. S2).

S3.3 Rank flux
In open ranking lists, elements flow into and out of the ranking (due to scores not being measured all the time), meaning that an element may not have a well-defined rank R t = 1, . . . , N 0 for all times t = 0, . . . , T − 1. We define rank flux F t as the probability that an element in (out of) the ranking list at time t − 1 leaves (enters) the ranking at time t, averaged over all elements in the ranking list 2 . For most open ranking lists, flux is roughly constant throughout time (upper rows in Fig. S3), i.e. F t does not considerably deviate from its mean over time, The Guardian readers (comm) Football scorers The constant trend in rank flux is sometimes disrupted by large deviations when a larger or smaller number of elements enter/leave the ranking list. In language datasets, F t has a decreasing trend over time (Fig. S3), arguably due to the exceptionally long observation period (centuries instead of days/years, see Table S1), or to a very heterogeneous sampling of words throughout observations. Closed ranking We define out-flux F − R as the probability that an element in rank R at time t − 1 leaves the ranking list at time t, averaged over all times t = 1, . . . , T − 1 (Fig. S4). Similarly, we define in-flux F + R as the probability that an element out of the ranking list at time t − 1 gets rank R at time t, averaged over time.

S3.4 Rank turnover
As mentioned in Section S1, N t is the number of distinct elements that have ever been in the ranking list up to time t (i.e. at any t = 0, . . . , t). Closed systems have N t = N 0 for any t, and open systems have N t > N 0 for some t > 0. The ranking list size N 0 is the initial condition of the time series N t , since in our initial observation we can only measure N 0 scores. Since N t counts the elements that have visited the ranking up until time t, it is a monotonically increasing function (N t ≤ N t ≤ N T −1 for t ≤ t). The value N T −1 is an observable proxy for the (unknown) size of the system that may increase with larger In open ranking lists, the time series o t increases from o 0 = 1 to o T −1 = N T −1 /N 0 (see Table S1), with In relatively open ranking lists, rank change is asymmetric: C is lower at the top of the ranking (low R) than at the bottom (high R). In progressively more closed ranking lists, rank change turns roughly symmetric in R, with low C also at the bottom of the ranking.

S3.5 Rank change
In order to measure the stability of rank dynamics, we introduce rank change C as the probability that elements in rank R at times t−1 and t are not the same, averaged over all times t = 1, . . . , T −1 (Fig. S6).
In C for both open and closed ranking lists.

S3.6 Rank inertia
We further analyze the flux of elements within the ranking list by dividing it into two regions via an arbitrary threshold c ∈ (0, 1): we define the top (+) of the ranking as all ranks lower than a threshold (R ≤ cN 0 ), and the bottom (−) as the rest of the elements in the ranking. We define the matrix element S ij t as the probability that an element in region i will move to region j after t observations (averaged over all compatible times t = 0, . . . , T − t − 1), with i, j ∈ {+, −}. The matrix element S ij t characterizes the flux of elements within the ranking list over time periods of size t. We focus on rank inertia S ++ t , the probability that an element stays in the top of the ranking list over a time t. Inertia S ++ t is a decreasing function of t for all datasets: while in some ranking lists the probability of staying in the top is roughly constant, for most systems it decays linearly, exponentially, or with a long tail (Fig. S7). In all

S4 Model of rank dynamics
Here we introduce and explore in detail a minimal model for the dynamics of ranking lists. We derive approximate master equations and their solutions for both the probability of an element changing rank and some of the rank measures introduced in Section S3. Our aim is to show how closely numerical simulations of the model (and their analytical approximations) emulate microscopic and macroscopic properties of the ranking dynamics of the datasets described in Section S2.

S4.1 Model definition
Our model is simple and corresponds to the null hypothesis that rank dynamics is driven by two processes: i) random displacements of elements across the ranking list, leading to Lévy flight-and diffusion-like movement in rank (tuned by a model parameter τ ∈ [0, 1]); and ii) random replacement of elements by new ones, leading to an increase in rank flux and turnover (tuned by a model parameter ν ∈ [0, 1]).
We assume that At each time step s we perform two independent rank updates. The first update, happening with probability τ , will cause Lévy flight-and diffusion-like dynamics. The second update, happening with probability ν, will lead to the replacement of elements in the ranking.
For the first update, we select an element uniformly at random, remove it from the list and place it in one of the N spaces to the right or left of the remaining N − 1 elements. More precisely, we choose rank R old = 1, . . . , N uniformly at random and take its element temporarily out of the system (thus making the rank change R → R − 1 for elements previously in ranks R = R old + 1, . . . , N ). Then we re-introduce this element to the rank R new = 1, . . . , N chosen uniformly at random (making the rank change R → R + 1 for elements previously in ranks R = R new , . . . , N − 1). In this way, elements initially between (and including) ranks R old and R new will change rank (if R old = R new ), while elements outside this range will keep their rank 4 (Fig. S8).
In the second update, we choose a rank R rep = 1, . . . , N uniformly at random and replace its element with a new element from outside the system, leaving the rest of the ranks untouched. The old element is removed from the system and ceases to have any dynamics, while the new element participates in the rank updates of time step s + 1 and beyond (Fig. S8). Each empirical dataset allow us to fix N 0 and T in the model (see Table S1). In what follows we will

S4.2 Displacement probability
At each time step of the dynamics, with probability ν∆r the element in rank R is replaced by a new element, thus leaving the system. Alternatively (with probability 1 − ν∆r), the element in rank R is displaced to a new rank X in the system with probability τ , either by being randomly selected with the displacement update rule of the model (Lévy flight), or because of the rank change of another element (diffusion). Like with R, we can define the normalized rank x = X∆r = ∆r, 2∆r, . . . , p, . . . , 1. We explore the temporal evolution of the model by formulating a master equation for the displacement probability P x,t , the probability that an element in the normalized rank r = R∆r moves to the normalized rank x = X∆r after a time t = s∆r.
We first consider the displacement probability after a single time step of the dynamics takes place (s = 1), P r x ≡ P x,∆r . An element in rank r (that is not replaced) moves to rank x with probability τ either by performing a Lévy flight (with probability L) or due to the rank change of another element (with probability D r x ), or stays in place with probability 1 − τ . Then we have with δ r x a Kronecker delta. Note that for τ = 1 and ν = 0 (in a closed ranking list with maximum number of displacements in a given time interval), the single-step displacement probability reduces to P r x = L + D r x . The Lévy probability L is straightforward to calculate: the element in rank r gets picked up by the dynamics with probability ∆r, and is placed back in rank x with the same probability ∆r, since both processes are uniformly random in rank space. Thus, Eq. (S5) explains our motivation to use the term Lévy flight: a Lévy flight is a random walk where displacement length follows a heavy-tailed probability distribution; since Eq. (S5) does not depend on the normalized displacement d = x − r, the probability L can be thought of as a discrete power law with exponent 0 (and cut-off due to finite system size) [38,39]. In other words, our model implements maximally heavy-tailed Lévy flights where all displacements are equally likely, beyond diffusion. If we would extend the model to allow for a non-uniform sampling of elements and ranks in the mechanisms of displacement and replacement, we expect that this exponent would become tunable.
The diffusion probability D r x is slightly more involved. The element in rank r moves one step to the right (x = r + ∆r) if we choose an element to its right (with probability 1 − x + ∆r) and place it back to its left (with probability x − ∆r). Conversely, the element in rank r moves one step to the left (x = r − ∆r) if we choose an element to its left (with probability x) and place it back to its right (with probability 1 − x). The element at rank r can stay in place (x = r) only if another element is picked and placed back to its left [with probability (x − ∆r) 2 ] or to its right [with probability (1 − x) 2 ]. Joining these expressions we obtain Note that x D r x = 1 − ∆r. In other words, the probability loss from the diffusive process in Eq. (S6) goes to the uniform contribution made by the Lévy term in Eq. (S5), meaning that the continuous limit of Eq. (S6) does not exactly correspond to a diffusion equation. Moreover, x P r x = 1 − ν∆r < 1 (for ν > 0), which implies a global loss of probability reflecting the fact that we replace elements in the dynamics.
Moving forward in the dynamics (t = s∆r ≥ 0), we write a recurrence relation for the displacement probability P x,t . We observe that, for an element in rank r to move to rank x after s time steps, it needs to have moved to an arbitrary rank x = ∆r, . . . , 1 over the first s − 1 time steps, and then from x to x in the last time step, i.e., with a Kronecker delta as initial condition (P x,0 = δ r x ), since at the start of the dynamics, an element in rank r can only be in that rank (x = r). Note that Eq. (S7) for s = 1 recovers P x,∆r = P r x in Eq. (S4), as expected.
Let's focus for a moment on the cumulative P t = x P x,t , the probability that, after a time t, the element in rank r moves to any other rank in the system. Summing up over x in Eq. (S7) gives the recurrence relation P t = (1 − ν∆r)P t−∆r , with initial condition P 0 = 1 (due to P x,0 = δ r x ) and P ∆r = 1 − ν∆r. Solving the master equation gives The approximation in Eq. (S8) comes from using the time scale relation t = s∆r, the binomial theorem, and the power series of the exponential, and becomes exact in the limit N → ∞. In other words, the probability of an element staying in the system under a dynamics of replacement (ν > 0) decays exponentially with time. When there is only element displacement (ν = 0), P t = 1 and the ranking list is indeed closed. Consistently, Eq. (S8) is equal to the probability that an element has not been replaced at any previous time, which decays geometrically or exponentially with time (in the discrete and continuous cases, respectively).
Inserting Eqs. (S4)-(S6) into Eq. (S7) we obtain a master equation for the displacement probability after s time steps, We can actually measure the displacement probability P x,t for t = 1 directly from empirical data as the probability that the element at rank R at time t − 1 moves to rank X at time t , averaged over all times t = 1, . . . , T − 1 (Fig. S9). The use of an average over time is supported by the stationarity of the empirical rank dynamics as seen in the roughly constant flux time series of most datasets (see Fig. S3).
As we will see in Section S4.3 below, the empirical displacement probability resembles either a flat Lévy sea or a diffusion peak, respectively corresponding to late and early times of the rank dynamics in the model. By calculating P x,t in the model for s = N (i.e. between times t = 0 and t = 1), we see that the model qualitatively reproduces the empirical displacement probability for most datasets and any rank [see Eq. (S22)].
As seen in Eq. (S8), the probability P t = x P x,t of an element staying in the system for ν > 0 decays exponentially with time, i.e. P t is not conserved. In order to approximately solve Eq. (S9) with a diffusion ansatz (see Section S4.3), we can renormalize P x,t by introducing the probability distribution which is indeed conserved, since x Q x,t = 1 for all t. As before, the approximation holds for t = s∆r for t = 1 as a function of normalized rank X/N0 for various R values (colored lines), calculated as the probability that the element at rank R at time t − 1 moves to rank X at time t , averaged over all times t . Px,t is shown for data (continuous lines) and the model of Section S4 with Eq. (S22) for s = N (dashed lines) (for parameter fitting see Table S2 and Section S5). Datasets are ordered from most open (upper row) to closed (lower row) according to mean rank flux F [Eq. (S2)]. In some datasets Px,t is roughly flat, reminiscent of a uniform Lévy sea, while in others the displacement probability resembles a diffusion peak (see Section S4.3). To smooth the curves, empirical data is passed through a Savitzky-Golay smoothing filter (with polynomial order 2 and filter window length ranging from N0/200 to N0/10). and is exact for N → ∞. Substituting Eq. (S10) into Eq. (S9) gives (S11) Eq. (S11) is a simplified master equation for the displacement probability that effectively decouples the dynamics of displacement (τ ) and replacement (ν) in the model. We can first look for an approximate solution Q x,t of the displacement dynamics governed by Eq. (S11), and then consider the replacement dynamics explicitly with P x,t = Q x,t e −νt .

S4.3 Approximation for displacement probability
Rather than solving Eq. (S11) explicitly, we derive an approximate expression for the renormalized displacement probability Q x,t that becomes more accurate as N → ∞. Considering the shape of Eq. (S4) (i.e. for s = 1), we propose an ansatz of Q x,t for any t separated into Lévy flight and diffusive components, (S12) As we will see below, the probability L t of moving across rank space in Lévy flights is rank-and displacement-independent and grows slowly in time (a uniform Lévy sea), while the probability D x,t of diffusing in rank space due to the movement of other elements is approximately Gaussian-distributed in x (a diffusion peak that widens and decreases in height as time goes by) (Fig. S10).
First, we write a recurrence relation for L t by noting that after a time t = s∆r, an element initially in rank r can only move due to diffusion within the 2s + 1 ranks around r, so Q x,t > L t for |x − r| ≤ t.
Conversely, the only way for an element to move more than s ranks to each side of r is by Lévy flight, meaning Q x,t = L t for |x − r| > t. Rewriting Eq. (S11) for |x − r| > t we obtain a recurrence relation for L t with initial condition L 0 = 0 (since at the start of the dynamics no Lévy flights have yet occurred), which consistently gives L ∆r = τ L as in Eqs. (S4)-(S5). Eq. (S13) can be directly solved and gives an exact expression for the probability of moving across rank space in Lévy flights, with t = s∆r and the approximation improving as N → ∞. In other words, the Lévy sea in Eq. (S14) is a uniform probability (both in ranks r, x and displacement d = x − r) that increases asymptotically from τ ∆r 2 to ∆r as t → ∞.
By inserting Eqs. (S12)-(S13) into Eq. (S11), we find a master equation for the probability D x,t of diffusing in rank space due to the movement of other elements, with initial condition D x,0 = δ r x . Just like Eq. (S11), Eq. (S15) is difficult to solve exactly for any t. We can, however, write an expression for A t = x D x,t , the area under the diffusion peak [disregarding the Lévy sea; see Eq. (S12)]. First, we know that x Q x,t = 1 due to its definition in Eq. (S10). Then, by plugging Eq. (S14) into this normalization condition we obtain with t = s∆r and the approximation getting better as N → ∞. In other words, L t = ∆r(1 − A t ).
Intuitively, the area A t of the diffusion peak decreases exponentially as time goes by, 'leaking probability' into a Lévy sea that increases in height with t.

S4.3.1 Continuum limit for diffusion peak
We may find an explicit, approximate solution for the diffusion peak D x,t by analyzing the continuum limit of Eq. (S15) as N → ∞. As mentioned in Section S4.1, in order to match the time scales between model (s) and data (t) we take t = s∆r with ∆r = 1/N . We also use the ansatz Q x,t = L t + D x,t with the Lévy sea L t given by Eq. (S14). Since the diffusion peak leaks probability into the Lévy sea as time goes by [see Eq. (S16)], we further propose the ansatz where D(x, t) is an unknown probability density, x is taken as a continuous variable, and the approximation improves for larger N . Inserting Eq. (S17) into Eq. (S15) we obtain a master equation for D(x, t), Eq. (S18) is a diffusion-like equation with a quadratic, rank-dependent diffusion coefficient αx(1 − x) and α τ ∆r (for large but finite N ). Note that the probability density D(x, t) decreases in time as the Lévy sea increases, so Eq. (S18) is not a standard diffusion equation. If we would modify Eq. (S6) to force A t = 1 by, say, adding a term ∆r in D r r , we would obtain a diffusion equation in which the rank-dependent coefficient is inside the first rank derivative. The initial condition D x,0 = δ r x in Eq. (S15) leads to the initial condition D(x, 0) = δ(x − x 0 ) (with x 0 ≡ r ∈ [0, 1] the initial rank of the element) and to the Dirichlet boundary conditions D(0, t) = D(1, t) = 0. Notice that Eq. (S18) is the Wright-Fisher equation [40], a continuous model of genetic drift that has been extensively studied from a mathematical point of view (see [41][42][43] and references therein).
We can solve Eq. (S18) exactly by separation of variables. Proposing the ansatz D(x, t) = u(x)v(t) leads to v(t) = v(0)e −λt with λ a separation constant, while u(x) needs to satisfy the eigenvalue equation (S19) By using the Frobenius method and the boundary conditions of Eq. (S18), we determine the allowed values of λ implicitly and the associated eigenfunctions u(x) as infinite series, over which the initial condition D(x, 0) can be expanded to obtain the particular solution of Eq. (S18) we are looking for.
However, it may be more instructive to find a closed-form approximation for D(x, t) that lets us qualitatively understand the behaviour of the displacement probability in rank and time. In a small enough interval around the initial rank x 0 = r, we can approximate x(1 − x) x 0 (1 − x 0 ) in Eq. (S18), which leads to the standard diffusion equation with diffusion coefficient αx 0 (1 − x 0 ). Then, the fundamental , the Gaussian distribution with mean x 0 and time-dependent standard deviation σ t = 2αx 0 (1 − x 0 )t, which spreads symmetrically around x 0 as time goes by 5 . This allows us to write where the term e −τ t is necessary for the Gaussian approximation to comply with Eq. (S16).
Finally, we derive an approximate expression for the displacement probability P x,t = Q x,t e −νt that becomes more accurate as N → ∞. Inserting Eq. (S20) into Eq. (S17), and using the ansatz of Eq. (S12) alongside Eq. (S14) we obtain with t = s∆r, or writing everything explicitly, Eq. (S21) captures the approximate behavior of the displacement probability P x,t intuitively: a Gaussian diffusion peak G(r, σ 2 t ) widening in time and leaking probability to the uniform Lévy sea L t , all regulated by an exponential loss in probability due to new elements. Eq. (S22) is a good approximation for P x,t even for relatively low N = 10 2 , as we can see by comparing with numerical simulations of the model described in Section S4.1 for s = N (i.e. t = 1) (Fig. S10). Eq. (S22) also shows that the model captures the displacement probability of the empirical datasets (see Fig. S9).

S4.4 Approximation for rank flux
Beyond the microscopic description of the dynamics given by the displacement probability P x,t , we also explore the temporal evolution of the model by approximating the rank measures introduced in Section S3 with closed expressions. We start with the mean rank flux F , measured for data as the probability that any element in the ranking list at time t − 1 leaves the ranking at time t, averaged over all recorded (observable) elements in the ranking list and over time [see Section S3.3 and Eq. (S2)]. To find F in the model, we define the time-dependent flux F t as the probability that a given element in rank r ≤ p leaves any of these ranks after time t = s∆r (either by displacement or by replacement). Following Eq. (S10), flux is given by where the step size in the sum is ∆r, i.e. x = ∆r, 2∆r, . . . , p.
We now find a master equation for the partial cumulative Q [a,b],t ≡ b x=a Q x,t for arbitrary a, b = ∆r, . . . , 1 and b > a 6 . Using Eq. (S11), summing over x, and changing dummy indices, we obtain Eq. (S24) is not closed on the cumulative Q [a,b],t , since it depends directly on the renormalized displacement probability (last two terms on the right side of the equation). To close the master equation we make the approximations Q a−∆r,t Q a,t and Q b+∆r,t Q b,t for all t (accurate as long as the diffusion peak of P x,t is far enough from a or b), which leads to We go back to the particular case of the flux F t . With a = ∆r, b = p, and the initial condition Additionally, we define F = 0 for p = 0 for consistency. Eq. (S27) is an analytical approximation for the probability that any element in the ranking list at time t−1 leaves the ranking at time t, according to our model. Intuitively, flux F increases with τ and ν, since more displacements or replacements make it more likely that elements will leave the ranking list. Conversely, larger p makes the ranking list longer (relative to system size) and thus less likely to lose elements due to out-flux. Eq. (S27) is a good approximation of F for all parameters considered, as we can see by comparing with numerical simulations of the model described in Section S4.1 (Fig. S11). Eq. (S27) also shows that the model reproduces the empirical flux of most datasets (see Fig. S3). As described in Section S3.3, in the data we can also measure what part of a ranking list contributes most to the flow of elements out of it by calculating the out-flux F − R : the probability that the element in rank R at time t − 1 leaves the ranking list at time t, averaged over all observed times. Similarly, in the model we define out-flux F − r,t as the probability that the element in the normalized rank r leaves the ranking list after time t = s∆r. We can compare F − R in the data with F − r,t for t = 1 (i.e. s = N ) in the model. An element contributes to out-flux either because it is replaced with a new element, or because it moves out of the ranking list via the displacement dynamics. Thus, where the displacement probability P x,t is given explicitly by Eq. (S22) in the limit of large N . Approximating the sum in Eq. (S28) by an integral in the interval x ∈ [p, ∞], changing variables and integrating by parts, we obtain reproduces the empirical out-flux of several datasets (Fig. S4). Rank out-flux increases with both τ and ν, since the displacement/replacement dynamics takes elements out of the ranking list. F − r,t is also relatively constant across the ranking list, except for bottom ranks (r ∼ p) where out-flux increases. Simulations are averages over 10 5 realizations (with T = 10) in a system of size N = 10 2 .

S4.5 Approximation for rank turnover
Here we derive an approximate, closed expression for the rank turnover o t = N t /N 0 , where N t is the number of distinct element that have been in the ranking list up to (and including) time t [see Eq. (S3)].
Similarly, in the model N t is the number of distinct elements that have been in the ranking list (r = ∆r, . . . , p) at any time step s ≤ s (with t = s∆r). Thus, our task is to find an explicit expression for N t . We start by introducing the probability p t that a randomly chosen element has been in the ranking list at any time t ≤ t, where M t is the number of distinct elements that have been in the whole system at any time t ≤ t.
Since the replacement dynamics adds one new element every time step with probability ν, on average we have M t N + νs, which leads to In order to find p t , we write a master equation for the probability p x,t that the element in any rank x of the model system 7 has been in the ranking list at any time t ≤ t, which we average to obtain We first note that the initial condition for p x,t is a step function over rank space. Within the ranking list, p x,t is always 1 and does not change in time, 7 x = ∆r, . . . , 1.
i.e., p x,t = p x,t−∆r = 1 for x = ∆r, . . . , p and t ≥ ∆r. In the rest of the system, p x,t increases from its initial condition 0. We now write the corresponding equation for elements outside of the ranking list, in a similar way to Eq. (S7). For an element in rank x to have ever been in the ranking list, it must have belonged to the ranking list at some point in its dynamics, and finally moved from some rank r to rank x in one time step. We can thus write With Eqs. (S4)-(S6) and some simplification we obtain a master equation for p x,t , where we have made the approximation 1 = p p,t p p+∆r,t in order to close the equation. We explicitly solve Eq. (S35) by recursion from the initial condition p via a geometric series, which leads to the asymptotic value Eq. (S37) clarifies the range of accuracy of the approximation p p+∆r,t 1. Depending only on the parameters of the model and as long as t is sufficiently large, if the fraction in the right hand side of Eq. (S37) is close to 1/p, then Eq. (S35) is an accurate master equation for p t . This typically happens for small τ or ν (see Fig. S13), or for p close to 1.
With Eq. (S31), Eq. (S3), and the explicit expression for p t in Eq. (S36), we finally obtain model that visit its ranking list as time goes by. Eq. (S38) recovers turnover in numerical simulations of the model (Fig. S13). In agreement with empirical measurements of o t (see Fig. S5), Eq. (S38) starts off with a concave shape from its initial condition o 0 = 1, and has a linear behavior for long times, where the slope of the line is regulated by the same fraction we see in Eq. (S37). Moreover, the mean turnover rate after t observations,ȯ t = (o t − o 0 )/t, has an asymptotic valueȯ ≡ȯ T −1 given bẏ

S4.6 Approximation for rank change
The approximation for the displacement probability P x,t in Eq. (S22) leads us directly to an explicit expression for the rank change C, defined for data as the probability that elements in rank R at times t − 1 and t are not the same, averaged over all times (see Section S3.5 and Fig. S6). In the model we define C t as the probability that an element located in rank r changes place after time t (we will finally compare C with C t for t = 1, i.e. s = N ). Since C t = 1 − P r,t , we have for t = s∆r and sufficiently large N . Eq. (S41) is a good approximation for rank change in the model that fails slightly at the extremes of the ranking list (Fig. S14). In agreement with empirical data (see is open (p < 1). Thus, in closed systems, the probability of changing rank is larger for elements in the middle of the ranking list, and lower at its extremes. In open systems, the bottom of the ranking list (r ∼ p) is also unstable.

S4.7 Approximation for rank inertia
Here we derive an approximate, closed expression for the inertia measure S ++ t introduced in Section S3.6.
Similarly to the case of empirical data, in the model we define the matrix element S ij t as the probability that an element in region i will move to region j after a time t = s∆r. Regions i, j are divided by an arbitrary threshold c = ∆r/p, . . . , 1 into the top (+) of the ranking (r = ∆r, . . . , cp) and the bottom (−) of the ranking (r = cp + ∆r, . . . , p). Inertia S ++ Intuitively, the probability of staying in the top of the ranking list decreases linearly with the flux of elements out of (or into) the ranking list. For p = 0 (a closed system) we get a slope β = 1, while for p = 1 and large τ (an open system with displacement dynamics) we get a slope regulated by the inertia threshold parameter, β c. This defines a region in (F, S ++ )-space where the model can reproduce pairs of (F, S ++ ) values coming from the empirical datasets, an area that becomes broader as c decreases.

S5 Fitting data with model
Here we describe the process of fitting the rank measures of all considered datasets with the minimal model of rank dynamics introduced in Section S4. First, we set N 0 , T and the time series N t with the empirical values of each dataset listed in Table S1. As described in Section S1 and Section S4. , which allows us to find τ * . The solution (τ * , ν * ) collapses to the universal curve τrνr = 1 when rescaling by the empirical values of ranking list size p, flux F , and mean turnover rateȯ. Data corresponds to the Scientists dataset [4] (see Table S1 and Section S2), but the graphical solution is qualitatively the same for all open ranking lists.
non-zero flux and out-flux, and turnover that increases in time. We find five strictly closed systems These relations lead to a non-linear system of equations for τ and ν, determined by the specific values of ranking list size, flux, and turnover of each dataset. Solving for ν (or τ ) leads to a transcendental equation in ν (or τ ) with no explicit solution. Still, we can gain insight by applying a graphical method and find the solution (τ * , ν * ) numerically (Fig. S16). From Eq. (S45a) we first note that, in order to ensure τ > 0, ν has to be in the interval pȯ < ν <ȯ. Sinceȯ is relatively low for all studied datasets (less than ∼ 0.3; see Table S2), we only need to consider small ν values in the fitting process. Within this range we find a single solution ν * by graphically equating the two sides of Eq. (S45b) (left panel in Fig. S16). This value can be inserted into Eq. (S45a) to find the remaining solution τ * (right panel in Fig. S16).
Out-flux F − R is relatively constant for most R values, and in open systems has a sharp increase for large R. This is only a good fit for relatively closed ranking lists, since in very open systems out-flux is a more gradually increasing function of R (Fig. S4). Turnover o t is concave and monotonically increasing in time, which follows empirical data (Fig. S5).

S5.1 Dynamical regimes of open ranking lists
Having explicit expressions in Eq. (S45) lets us further analyze the solution (τ * , ν * ). We introduce the rescaled parameters which allows us to rewrite Eq. (S45a) for sufficiently small ν (relative to p and 1 − p), τ r ν r 1.
Eq. (S48) shows that, when fitting the model to datasets with given values of p, F , andȯ, the probabilities τ and ν can be rescaled to collapse into a universal curve, a power law with exponent 1 (Fig. S16).
In order to reproduce the dynamics of empirical ranking lists, is it enough to consider an inverse relationship between the rates of displacement and replacement of elements in the model, meaning that the changes rank by diffusion  or is replaced, where W levy + W diff + W repl = 1.
In the most open systems we study, elements mostly change rank via long jumps, forming a Lévy walk regime where W levy > W diff , W repl (Fig. S17). Here, long-range rank changes take elements in and out of a small ranking list within a big system (low p), thus generating large mean flux F . The rest of the datasets belong to a diffusion regime with W diff > W levy , W repl . A local, diffusive rank dynamics is the result of elements smoothly changing their scores and overcoming their neighbors in rank space. Although not observed in empirical data, the model predicts a third regime dominated by replacement (W repl > W levy , W diff ). The curves in Fig. S17 show how close a ranking list might be to a change of regime. If a dataset has values of p andȯ that make it fall close to the boundaries of these regimes, small changes in the parameters regulating the displacement (τ ) and replacement (ν) can change the dynamical regime of the ranking list.

S5.2 Fluctuations, estimation bias, and goodness of fit
The relationship between rescaled parameters τ r and ν r in Eq. recover the empirical values of system-wide quantities (F ,ȯ, and inertia S ++ ), showcasing the ability of the model to reproduce rank dynamics in data. We further quantify the goodness of fit by showing that, for most systems considered, data is closer to the universal curve [Eq. (S48)] than model simulations.
We start by taking a dataset with given values of T , p = N 0 /N T −1 , and fitted parameters τ and ν (see Tables S1-S2). We first run a large number of realizations of the model with the same values of N 0 , T , τ , and ν, plus an arbitrary but smaller system size N sim < N T −1 , leading to a certain number of elements ever seen in the ranking list of each simulation, which we denote by N T −1,sim . Since the model We compare the distributions of p sim , τ sim , and ν sim over model simulations with the fitted parameters p, τ , and ν of open ranking lists by calculating the differences ∆p = p sim − p, ∆τ = τ sim − τ , and ∆ν = ν sim − ν (Fig. S18). We obtain ∆p ∆ν 0 for all datasets (apart from hyenas, where ν is slightly underestimated), meaning that the fitting process consistently recovers the values of p and ν in simulations, without bias. Fluctuations around ∆p , ∆τ , and ∆ν (measured by the spread of their associated distributions) are only noticeable for some of the smallest datasets (Enron emails, hyenas, Nascar drivers, cities [RU], universities, and countries). Even when accounting for fluctuations, we find a small bias in τ (| ∆τ | < 0.1) for several datasets [The Guardian readers (comm), hyenas, languages, scientists, cities (RU), universities, countries, and chess], indicating that the fitting process systematically under-(∆τ < 0) or over-(∆τ > 0) estimates the displacement probability in simulations, particularly in systems with low N T −1 or T .
Despite the small bias in τ (Fig. S18), numerical simulations with fitted model parameters recover aggregate properties of empirical data, which we quantify via the differences ∆F = F sim −F , ∆ȯ =ȯ sim − o, and ∆S ++ = S ++ sim − S ++ (Fig. S19) In order to further quantify the goodness of fit, we explore how the rescaled fitted parameters and p,ȯ, and fit τ and ν, the rescaled parameters τ r and ν r in Eq. (S47) have the same relative difference from the expected universal behavior (τ r ν r = 1) alongside either axis, τ r − 1/ν r 1/ν r = ν r − 1/τ r 1/τ r = τ r ν r − 1.
We measure τ r ν r −1 for all datasets, as well as the distribution of the corresponding quantity τ sim r ν sim r −1 over model simulations, calculated with p sim ,ȯ sim , and the bootstrapped fit τ sim and ν sim 8 (Fig. S20).
We do not expect rescaled parameters to follow the universal curve perfectly, due to the approximations used to derive F andȯ in the model [Eq. (S27) and Eq. (S40); see Sections S4.4-S4.5], as well as the assumption of small ν in the derivation of Eq. (S48) itself (see Section S5.1 and Fig. S16). Still, in all datasets we find good agreement for fitted model parameters (|τ r ν r − 1| < 0.5) and bootstrapped estimates (| τ sim r ν sim r − 1| < 0.5). For most systems, the dataset has τ r and ν r closer to the universal curve than most simulations of the model itself, suggesting than empirical ranking lists follow Eq. (S48).
Exceptions are The Guardian readers (recc), Enron emails, Nascar drivers (Busch), and hyenas. We consider that these systems do not follow the proposed universal behavior either because of small sample size (leading to data noise), or due to a systematic failure of the model to recover their dynamics.
Finally, it's worth noting that it would be ideal to have a more formal statistical approach to model and (ii) the probability of element displacement in Eq. (S21). The maximum likelihood equations (for τ and ν) coming from maximizing the log-likelihoods of these two distributions are trascendental and cannot be solved analytically. We can, however, solve the equations numerically and compare their solutions to the fitted τ and ν from our 'mean-field' fitting process. For most datasets, mean-field parameter estimates lead to lower mean-squared-error deviations between the fitted model and several measures in the data (rank flux, turnover, out-flux, change probability, rank diversity, and inertia). This means that our fitting process is able to recover aggregate features of empirical ranking data, despite it being based on asymptotic approximations of mean behavior and showing a small bias in estimating τ . As added value, the mean-field fitting process is more analytically tractable, computationally faster to obtain, and provides insight on the balance between amounts of element displacement and replacement seen in data.

S5.3 Effect of subsampling
Our model of ranking dynamics suggests that empirical ranking lists belong to one of several dynamical regimes constrained by the universal curve in Eq. (S48). Here we explore the role of sampling rate (related to the time between observations of the ranking, see in Table S1) in the regime occupied by a system. We show that subsampling observations typically makes systems go downwards along the universal curve (right panel in Fig. S16 and Fig. 3 in main text), moving, for example, from a regime with many Lévy walks to one more driven by diffusion.
For a synthetic or empirical system with T observations, we implement subsampling by only considering the ranking list at times t = 0, k, 2k, . . . , T eff − 1 and discarding the rest of the data, with k = 1, 2, . . . , T − 1 a tunable parameter. The subsampling process decreases the number of observations to T eff = T /k ≤ T , such that k = 1 recovers the original dataset (T eff = T ) and k = T − 1 corresponds to the most subsampled system possible (T eff = 2). We then calculate the (k-dependent) subsampled rank flux F k , turnover rateȯ k , and inertia S ++ k , and compare them with their values F ,ȯ, and S ++ in the original ranking list. Finally, we fit the model to the subsampled data [via Eq. (S45) or Eq. (S46)] for each k value, obtaining the subsampled relative list size p k , displacement τ * k and replacement ν * k , in contrast with the original fit p, τ * , and ν * .
We can gain some insight on the role of subsampling by first analyzing synthetic ranking lists coming from the model itself. Consider a model dynamics (with parameters N , τ , and ν) and its subsampled version (with effective parameters τ k and ν k for given k) (Fig. S21). The probability that an element is not replaced between two subsampled observations of the original dynamics, (1 − ν∆r) kN , must be equal to the probability of not being replaced between consecutive observations of the subsampled dynamics, (1 − ν k ∆r) N . A similar argument applies to τ and τ k . Taking the limit N → ∞ and the first-order power series of the exponential, we obtain meaning that subsampling increases the rates of displacement and replacement in the model, almost linearly for small enough τ and ν. Eq. (S53) allows us to write approximate analytical expressions for the (for threshold c = 0.5), as well as fitted parameters p k , τ * k and ν * k as a function of subsampling parameter k, relative to their original values (F ,ȯ, S ++ , p, τ * , and ν * ) for fixed τ , ν, and p. Results are both numerical simulations of the model (dots) and their corresponding analytical approximations (lines), with τk and νk given by Eq. (S53). Subsampling increases the rates of displacement and replacement, leading to larger flux and turnover, and lower inertia. We only show k values such that Teff /T > 0.05. Simulations are averages over 10 realizations (with T = 10 3 ) in a system of size N = 10 3 .
subsampled rank flux F k , turnover rateȯ k , and inertia S ++ k by replacing (τ, ν) with (τ k , ν k ) in Eq. (S27), Eq. (S38), and Eq. (S43), respectively. As k grows, larger rates of displacement and replacement lead to an increase in flux and turnover, and a corresponding lower probability to stay in the top of the ranking.
Since some elements enter and leave the ranking list between subsampled observations, the number of elements ever seen in the list tends to decrease, leading to a slightly larger p k (Fig. S21). We note that the approximation of Eq. (S53) might become less accurate for some parameter values and as k increases.
We perform the same subsampling process for all open empirical ranking lists, and further calculate the rescaled parameters τ r and ν r for varying k [by using Eq. (S47) with the subsampled quantities τ * k ν * k , p k , andȯ k instead of their original values] (Fig. S22). As subsampling increases, most datasets go downwards along the universal curve of Eq. (S48) in (ν r , τ r )-space, moving from a regime with a certain amount of Lévy walks to one more driven by diffusion. Some datasets do not follow this trend [Nascar drivers, hyenas, cities [RU], companies, universities], which follows from the condition ν * k ∼ 0 in Eq. (S48) not being fulfilled for large enough k (see Fig. S21). We also notice some overlap between these datasets and the ones showing both large fluctuations and distance to the universal curve in Section S5.2 ( Fig. S18 and Fig. S20), further suggesting that only a few empirical ranking lists have dynamics not captured by the model and its universal curve.

S6 On the heterogeneity of rank dynamics
Our results suggest that rank dynamics are heterogeneous: some elements change their rank slower than others, even when the dynamics of the elements themselves (i.e. the scores) is homogeneous. Soccer  So why are rank dynamics heterogeneous? A possibility is that systems require additional mechanisms for maintaining a homogeneous temporality. Still, are there advantages to such a heterogeneity? When considering open systems in close interaction with their surroundings (and for which the entirety of the associated ranking has a physical meaning, e.g. in the case of a genetic code), there is evolutionary advantage in preserving the functionality of the most essential elements (to maintain robustness [44]), while allowing for fast variability of less crucial components (conferring adaptivity [45]). An understanding of the balance between stability and variability has so far been mostly constrained to critical phase transitions in homogeneous models [46][47][48]. We pose the following hypothesis: if we assume varying timescales of dynamical behavior according to rank (with elements at the top of a ranking list changing more slowly than the rest), systems would benefit by having robustness and adaptability at the same time, independently of critical parameters: "slower" elements would provide robustness, while "faster" elements would provide adaptivity. The ability of a simple model to reproduce rank dynamics in a wide variety of phenomena, regardless of their domain, suggests that such balanced dynamics can be achieved with random rank change. Thus, complex systems may need much less to be evolvable than previously thought (on the basis of homogeneous timescales only), implying that random variation of natural selection, if heterogeneous, is enough to produce the complex adaptations seen in evolutionary biology and computer science [49]. This hypothesis remains to be explored further, but it might be useful to guide the interpretation of the results presented here.