{'fraction_non_alphanumeric': 0.03724237010295987, 'fraction_numerical': 0.01256147201296668, 'mean_word_length': 4.533992457445724, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 9, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Emergence of Hierarchy in Networked Endorsement Dynamics

Many social and biological systems are characterized by enduring hierarchies, including those organized around prestige in academia, dominance in animal groups, and desirability in online dating. Despite their ubiquity, the general mechanisms that explain the creation and endurance of such hierarchies are not well understood. We introduce a generative model for the dynamics of hierarchies using time-varying networks in which new links are formed based on the preferences of nodes in the current network and old links are forgotten over time. The model produces a range of hierarchical structures, ranging from egalitarianism to bistable hierarchies, and we derive critical points that separate these regimes in the limit of long system memory. Importantly, our model supports statistical inference, allowing for a principled comparison of generative mechanisms using data. We apply the model to study hierarchical structures in empirical data on hiring patterns among mathematicians, dominance relations among parakeets, and friendships among members of a fraternity, observing several persistent patterns as well as interpretable differences in the generative mechanisms favored by each. Our work contributes to the growing literature on statistically grounded models of time-varying networks.

Many social and biological systems are characterized by enduring hierarchies, including those organized around prestige in academia, dominance in animal groups, and desirability in online dating. Despite their ubiquity, the general mechanisms that explain the creation and endurance of such hierarchies are not well understood. We introduce a generative model for the dynamics of hierarchies using time-varying networks in which new links are formed based on the preferences of nodes in the current network and old links are forgotten over time. The model produces a range of hierarchical structures, ranging from egalitarianism to bistable hierarchies, and we derive critical points that separate these regimes in the limit of long system memory. Importantly, our model supports statistical inference, allowing for a principled comparison of generative mechanisms using data. We apply the model to study hierarchical structures in empirical data on hiring patterns among mathematicians, dominance relations among parakeets, and friendships among members of a fraternity, observing several persistent patterns as well as interpretable differences in the generative mechanisms favored by each. Our work contributes to the growing literature on statistically grounded models of time-varying networks.
introduction Hierarchies-stable sets of dominance relationships among individuals [Fushing et al., 2011;Hobson et al., 2018;Hobson and DeDeo, 2015]-structure many human and animal societies. Among animals, hiearchical rank may determine access to resources such as food, grooming, and reproduction [Holekamp and Strauss, 2016]. Among humans, rank shapes the epistemic capital and employment prospects of researchers [Clauset et al., 2015;Morgan et al., 2018], susceptibility of adolescents to bullying [Garandeau et al., 2014], messaging patterns in online dating [Bruch and Newman, 2018], and influence in group decision-making [Cheng and Tracy, 2014].
A central question concerns how enduring hierarchies shape and are shaped by interactions between individuals. Empirical studies have indicated the presence of a winner effect: an individual who participates in a favorable interaction, such as winning a fight or receiving an endorsement, increases their likelihood of being favored in future interactions [Chase et al., 1994;Hogeweg and Hesper, 1983]. Both theoretical work [Bonabeau et al., 1995[Bonabeau et al., , 1996Hemelrijk, 1999;Ben-Naim and Redner, 2005;Miyaguchi et al., 2020;Pósfai and D'Souza, 2018;Hickey and Davidsen, 2019;Vehrencamp, 1983;Sánchez-Tójar et al., 2018] and controlled experiments in humans [Salganik et al., 2006] suggest that winner effects are sufficient (though not necessary) to form stable hierarchies. Mechanistic explanations of winner effects vary. A common approach postulates that each individual possesses an intrinsic strength, which may depend on factors such as size, skill, or aggression levels. For instance, physiological mechanisms such as changes in hormone levels following confrontational interactions [Mehta and Prasad, 2015] can alter an individual's strength, causing the strong to get stronger.
However, intrinsic strengths are not necessary to produce winner effects. If a politician endorses a rival candidate, the latter does not become intrinsically more fit for office; instead, the endorsee builds support for their candidacy that may lead to future endorsements. The fame of the endorser is key: the better-known the endorser, the more valuable the endorsement. We refer to such prestige by proxy as transitive prestige. Since transitive prestige enables hierarchical rank to flow through interactions between individuals, networks provide a natural lens through which to study its role. Recent empirical studies have emphasized the networked nature of hierarchy in biological and social groups [Shizuka and McDonald, 2015;Ball and Newman, 2013;Pinter-Wollman et al., 2014;Hobson and DeDeo, 2015;Hobson et al., 2018]. Several theoretical studies [König and Tessone, 2011;König et al., 2014;Bardoscia et al., 2013;Krause et al., 2013] have also investigated reinforcing hierarchy using time-varying network models called adaptive networks [Sayama et al., 2013;Porter, 2020]. In this class of models, edges, representing interactions, evolve in response to node states and vice versa. Edges tend to accrue to important or highly central nodes, leading to self-reinforcing hierarchical network structures. Despite their recent uses, adaptive networks are often difficult to analyze analytically or compare to empirical data.
We present a novel and flexible adaptive network model of social hierarchy that addresses these challenges. Winner effects in our model are driven entirely by social reinforcement rather than intrinsic strengths. We allow arbitrary matrix functions to determine rank or prestige of nodes in the network, and introduce parameters governing the behavior of individuals in response to rank. A key feature of our model is that it is amenable both to mathematical analysis and to statistical inference. We analytically characterize a critical transition separating egalitarian and hierarchical model states for several choices of ranking function. We also explore hierarchical patterns in four biological and social data sets, using our model to perform principled selection between competing ranking methods in each data set and highlight persistent macroscopic patterns. We conclude with a discussion of potential model extensions and connections to recent work on centrality in temporal networks.

modeling emergent hierarchy
In our adaptive network model, new directed edges are formed based on existing, node-based hierarchy, after which they decay over time. We conceptualize a directed edge i → j as an endorsement, in which i affirms that j is fit, prestigious, or otherwise of high quality. For example, endorsements could capture contests won by j over i, retweets of j by i, or comparisons in which a third party ranks j above i. We collect endorsements in a weighted directed network on n nodes summarized by its adjacency matrix A ∈ R n×n , where entry a ij is the weighted number of interactions i → j. The matrix A evolves in discrete time via the iteration Here, the update matrix ∆(t) contains new endorsements at time t. The memory parameter λ ∈ [0, 1] represents the rate with which memories of old endorsements decay; the smaller the value of λ, the more quickly previous endorsements are forgotten.
The new endorsements in ∆(t) depend on previous endorsements through a ranking of the n nodes, which we call the score vector (or simply score) s ∈ R n . The score vector is the output of a score function σ : A → s ∈ R n , which may be any rule that assigns a real number to each node.
We consider three score functions chosen for analytical tractability and relevance in applications. Let D in and D out be diagonal matrices whose entries are the weighted in-and out-degrees of the network, i.e., D in ii = j A ij and D out ii = j A ji . First, the Root-Degree score is the square root of the in-degreethe weighted number endorsements-of each node i, defined as s i = D in ii . The Root-Degree score function does not model transitive prestige, since only the number of endorsements is considered, not the prestige of the agents from which they come. Second, PageRank [Brin and Page, 1998] is a recursive notion of rank in which high-rank nodes are those whose endorsers are numerous, and themselves high rank. The foundational algorithm used by Google in ranking webpages, PageRank computes a value for each node interpretable as the proportion of time that a random surfer following the network of endorsements would spend on that node. We define PageRank score s as the PageRank vector of A T , which is the unique solution to the system up to scalar multiplication. Here, α p ∈ [0, 1] is the so-called teleportation parameter, for which we use the customary value α p = 0.85. We normalize the PageRank vector so that e T s = n, where e is the vector of ones. Finally, SpringRank [De Bacco et al., 2018] is another recursive definition of rank in which endorsers are ranked one unit below endorsees, and disagreements are resolved using an analogy to a physical system of springs: the ranking of nodes minimizes the total energy of the system. Mathematically, the SpringRank score s is the unique solution to the linear system [De Bacco et al., 2018] with the identity matrix I and a regularization parameter α s > 0 which ensures the uniqueness of s. Unlike the Root-Degree score, both PageRank and SpringRank scores model transitive prestige: the impact of an endorsement depends on the prestige of the endorser. These three score functions can all be interpreted as rankings or centrality measures, although this property is not required of score functions in our model. Given score vector s, new endorsements ∆ are chosen using a random utility model, a standard framework in discrete choice theory which has recently been applied in models of growing networks [Overgoor et al., 2019]. At time step t, node i is selected uniformly at random. We suppose that endorsing j has utility u ij (s) for i, which depends on the current scores. In this work, we focus on utilities of the functional form where we generally assume that β 1 > 0 and β 2 < 0. The parameter β 1 captures a preference for prestige; a positive value of β 1 indicates a tendency to endorse others with high scores. The parameter β 2 captures a preference for proximity; a negative value of β 2 indicates a tendency to endorse others with scores relatively similar to their own. Many other choices of utility functions are possible; we prove a stability theorem for a large class of these functions in Appendix A.
In the random utility model, node i observes all possible utilities subject to noise. Traditionally, this noise is chosen to be Gumbel-distributed, in which case the probability that endorsing j yields the greatest utility is given by the multinomial logit [Train, 2009] We collect m ∈ N endorsements in an update matrix ∆, where the entry ∆ ij gives the number of times that i endorses j in the time step. More complex random utility models can lead to more realistic structures in networks with a growing number of nodes [Gupta and Porter, 2020]; we do not pursue these complications here because our model does not focus on network growth, and because these complications obstruct analytical insight. Equation (5) can also be derived from an alternative model in which node i makes a randomized choice among n nodes to endorse. In this model, the option to endorse j is assigned a deterministically-observed weight proportional to e uij (s) . In this case, β 1 and β 2 signify inverse temperatures that tune the degree FIGURE 1 -Schematic illustration of our model dynamics. Nodes are initialized at time t = 1 with a set of pre-existing endorsements logged in A (solid arrows) and the score s = σ(A) is computed (vertical axis). Then, a new edge logged in ∆ is added (dashed line). In the next time step t = 2, old interactions decay by a factor of λ (grey arrows). The new endorsement and decay of previous endorsements lead to an updated score function, which then informs the next time step.
of randomness in this choice, with lower values corresponding to greater randomness. Although this alternative model-in which node i makes a noisy choice between deterministically-observed utilities-and the random utility model-in which node i makes a deterministic choice between noisily-observed utilities-are mathematically equivalent, the two formulations can lead to different interpretations of system behavior. In the case of institutional faculty hiring discussed below (see Hierarchies in Data), the random utility model assumes that a hiring committee makes imperfect observations of the utilities of the institutions from which they could hire, and then deterministically chooses the highest of these imperfectly-observed qualities. In contrast, the alternative framework assumes that the committee makes a perfect observation of the utilities, but then chooses among them with some degree of randomness, which may reflect dissension on the hiring committee, search-specific priorities, among other factors. Equations (1) and (5) capture key features of our model. First, the dynamics in Eq. (1) imply that past interactions decay geometrically at rate λ. This global, gradual decay contrasts with another rank-based relinking model in which single edges fully disappear within each time step [König and Tessone, 2011]. Second, Eq. (5) implies that the likelihood of a node being endorsed at a given time step depends only on the distribution of previous endorsements and not on intrinsic strength or desirability. Those who receive more endorsements and therefore obtain higher scores are more likely to be endorsed in the future-a mechanistic instantiation of winner effects via social reinforcement. Figure 1 schematically illustrates model dynamics with m = 1 endorsement per time step. At time t = 1, the model is initialized with a small number of endorsements logged in A. The score function takes A as an input and outputs the score vector s, which in turn determines a new interaction according to Eq. (5). Logged in ∆, this new interaction is weighted by 1 − λ and added to the previous endorsements, which are discounted by λ. This process repeats over time with new endorsements gradually replacing old ones in the system's memory, sequentially updating the score vector s. Figure 1 also depicts in stylized fashion the operation of both a winner effect (β 1 > 0), in which endorsements tend to flow in the direction of increasing score, and a proximity effect (β 2 < 0), in which endorsements tend to flow between nodes of similar scores. The net effect is that most endorsements are "short hops" up the hierarchy. As we will discuss, this is a common pattern in empirical data.
Despite its simplicity, the model displays a wide range of behaviors. To observe them, we define a rank vector γ, whose jth entry γ j = n −1 i p ij gives the likelihood that a new endorsement flows to j. We say that the system state is egalitarian when all ranks γ j are equal and hierarchical otherwise. Figure 2 illustrates representative behaviors when the SpringRank score is used. When β 1 is relatively small, winner effects are overtaken by noise, and the system settles into an approximately egalitarian state ( Fig. 2a,b). When β 1 is relatively large, persistent hierarchies emerge ( Fig.2c-f). Moreover, the distribution and stability of ranks depend on the strength of proximity effects, modeled by the quadratic term in the utilities. For β 2 = 0 (no proximity preference), a single node garners more than half of endorsements in a hierarchy with significant fluctuations (Fig. 2c,d). Adding a proximity preference leads to a marginally more equitable hierarchy with ranks that are nearly constant in time (Fig. 2e,f).

the long-memory limit
The behavior observed in Fig. 2 suggests the presence of qualitatively distinct regimes depending on prestige preference β 1 . For small β 1 (Fig. 2a), the winner effect is weak and approximate egalitarianism prevails. For larger β 1 , a stronger winner effect enforces a stable hierarchy. We characterize the boundary between these regimes analytically in the long-memory limit λ → 1 by defining a function f which is analogous to a deterministic time-derivative for the dynamics of our discrete-time stochastic process. Let where the expectation is taken with respect to ∆. If f (s, A) = 0 for all A, the score vector s is a fixed point of the model dynamics in expectation. Our choices of Root-Degree, PageRank, and SpringRank score functions admit closed-form expressions for f , allowing us to analytically derive the conditions for the stability of egalitarianism in the limit of long memory.
Theorem 1. For each of the Root-Degree, PageRank, and SpringRank score functions, f has a unique egalitarian root. This root is linearly stable if and only if β 1 < β c 1 , where In Appendix A, we prove Theorem 1, as well as a generalization to arbitrary smooth utility functions. In each case, the proof of uniqueness exploits the FIGURE 3 -Bifurcations in models with Root-Degree, PageRank, and SpringRank score functions with β 2 = 0 and m = 1 update per time step. Points give the value of the rank vector γ averaged over the final 500 time-steps of a 5 × 10 4 -step simulation with n = 8 nodes, memory parameter λ = 0.9995, and varying β 1 specified by the horizontal axis. Solid curves show stationary points of the longmemory dynamics obtained by numerically solving the equation f (s, A) = 0, subject to the restriction that nodes separate into two groups with identical ranks in each. Black curves are linearly stable, while grey curves are unstable. Stability was determined by studying the spectrum of the Jacobian matrix of f . Vertical lines give the critical value β c 1 at which the egalitarian solution becomes linearly unstable according to Theorem 1. Parameters: αp = 0.85, αs = 10 −8 .
algebraic structure of the score function, and the critical value β c 1 is obtained via the linearization of f about the egalitarian state. Interestingly, only β 1 plays a role in the stability of the egalitarian root. While proximity preference β 2 does not determine where the hierarchical regime begins, it does influence the structure of and the transient dynamics toward nonegalitarian equilibria ( Fig. 2(e),(g)). Figure 3 illustrates the destabilization of egalitarianism predicted by Theorem 1 in the case of n = 8 nodes. Although not required by Theorem 1, we fix β 2 = 0 for simplicity. Curves show fixed points of the model dynamics in the long-memory limit. We show only fixed points in which nodes separate into two groups, each of which have identical rank. For β 1 < β c 1 , the egalitarian regime is stable and the long-run state deviates from egalitarianism only slightly. For β 1 > β c 1 , in contrast, the long-run state switches to an inegalitarian, stable fixed point.
In the Root-Degree and PageRank models, there is a single stable inegalitarian equilibrium with one node absorbing nearly all endorsements ( Fig. 3a,b). Interestingly, there is a bistable regime in which both egalitarian and inegalitarian states are attracting. Whether the system converges to one or the other depends on initial conditions. The SpringRank model displays qualitatively distinct behavior (Fig. 3c,d). Past β c 1 , we observe staggered multistable regimes. As β 1 increases, equilibria with multiple elite (i.e., highly ranked) nodes become sequentially unstable until eventually only a single elite node remains. The long-term behavior of the system again depends on initial conditions, but now there are many more possible stable states. This behavior would seem to make the SpringRank score function especially appropriate for modeling empirical systems with multiple distinct hierarchical regimes and sensitivity to initial conditions, an intuition which we confirm empirically in the following section.

hierarchies in data
In addition to being amenable to analytical treatment, our model has a tractable likelihood function, described in Appendix B. This allows us to study hierarchical structures in empirical data using principled statistical inference. The likelihood function not only supports maximum-likelihood parameter estimates of λ, β 1 , and β 2 but also enables direct comparisons of different score functions in a statistically rigorous framework: score functions with higher likelihoods provide more predictive low-dimensional summaries of observed interactions. This in turn allows us to explore the relative value of competing mechanistic explanations of observed data.
Several mathematical features of the model facilitate the exploration of real data. First, the predictive distribution Eq. (5) is in the linear exponential family, making the estimation of β a convex optimization with a unique solution. Second, the estimation problem inλ is in general nonconvex, but can be tractably solved via first-order optimization methods with multiple starting points. Finally, while model likelihoods evaluated on training data may in principle be inflated due to overfitting, our model uses only three parameters to fit hundreds or thousands of observations, suggesting that overfitting is not a major concern.
We conducted a comparative study of model behavior on four data sets: an academic exchange network in math, two networks of parakeet interactions, and a network of friendships among members of a fraternity. The Math PhD Exchange data set is extracted from The Mathematics Genealogy Project [North Dakota State University Department of Mathematics, 2020;Taylor et al., 2017]. Nodes are universities. An interaction i → j at time t occurs when a mathematician who received their degree from university j at time t supervises one or more PhD theses at university i. This event is a proxy for university i hiring a graduate from university j at a time near t. We view this as an endorsement by j that graduates of i are of high quality [Clauset et al., 2015]. We restricted our analysis to the activity of the 70 institutions that placed the most graduates between 1960 and 2000. Doing so helped to avoid singularities produced by institutions with no placements early in the time period and to minimize temporal boundary effects associated with the beginning and end of data collection.
The two Parakeet data sets [Hobson and DeDeo, 2015] record aggression events in two distinct groups of birds studied over four observation quarters (weeks). An interaction i → j at time t occurs when parakeet i loses a fight to parakeet j in period t. Since there are just four observation periods, estimates of the memory parameter λ should be approached with caution.
Lastly, the Newcomb Fraternity data set was collected by the authors of Nordlie [1958] and Newcomb [1961] and accessed via the KONECT network database [Kunegis, 2013]. The data set documents friendships among members of a fraternity at the University of Michigan. Each week during a fall semester, excluding a week for fall break, each of 17 cohabiting brothers ranked every other brother according to friendship preference, with ranks 1 and 16 referring to that brother's most and least preferred peers, respectively. An endorsement i → j is logged when brother i ranks j among his top k = 5 peers (small changes to k did not significantly alter the results). While friendship is often viewed as a symmetric relationship, expressed friendship preferences may be asymmetric [Carley and Krackhardt, 1996].
We studied these data using the Root-Degree, PageRank, and SpringRank score functions. Table 1 summarizes our results, including parameter estimates; standard errors (obtained by inverting the numerically-calculated Fisher information matrix); and optimized log-likelihoods for each combination of score and data set. Several features stand out. In all four data sets and across all three score functions, we findβ 1 > 0 andβ 2 < 0. This suggests a persistent pattern in time-dependent hierarchies: while endorsements do flow upward (β 1 > 0), nodes are more likely to endorse those close to them in rank (β 2 < 0). Endorsements tend to flow a few rungs up the ladder-not directly to the top. The reasons for this pattern likely vary across data set. In the Math PhD Exchange, this may indicate that low-ranked schools struggle to recruit graduates of high-ranked ones due to a limited supply of elite candidates. In parakeet populations, proximal aggression may facilitate inference of dominance hierarchies through transitive inference [Hobson and DeDeo, 2015]. In Newcomb's Fraternity, we postulate that implicit social norms may encourage friendships between those of similar standing. Similar results have been reported in static social network data among adolescents [Ball and Newman, 2013]. Thus, while we do not attribute this pattern  For each data set, the largest log-likelihood L is indicated in bold. All parameter estimates are statistically distinct from zero at 95% confidence. N gives the total number of interactions in the data. See Fig. 9 for simulated trajectories with the inferred parameters.
in the parameter estimates to a "universal" mechanism, we suggest its persistence as an interesting observation worthy of future study. Because different score functions capture distinct qualitative features of the data, quantitative comparisons yield insights into the generating mechanisms at work. In general, parameters from models using differing score functions should not be directly compared, since these parameters are sensitive to the scale of the score vector. However, we can compare models on the basis of their likelihoods. In the Math PhD Exchange, the Root-Degree model was strongly favored over either SpringRank or PageRank. In the context of this data set, the Root-Degree score is a measure of faculty production: a school that places more candidates has a higher score, regardless of the prestige of the institutions at which the candidates land. The strong fit from the Root-Degree score is consistent with previous findings that raw faculty production plays a major role in structuring the hierarchy of academic hiring within computer science, business, and history [Clauset et al., 2015]. As Clauset et al. [2015] note, transitive prestige also plays an important role. It would be of significant interest to extend our study to include multiple score functions, enabling an inferential analysis of the relative roles of production and transitive prestige.
In contrast, the SpringRank score was favored by large margins in both Parakeet data sets and by a smaller margin in the Fraternity data set, suggesting that transitive prestige plays a more prominent role. Among Parakeets, it may matter not only how many confrontations one wins but also against whom, with victories over high-ranking birds counting more towards one's own prestige. This finding is consistent with those of Hobson and DeDeo [2015], who found, using different methodology, that parakeet behavior suggests the ability to draw sophisticated, transitive inferences about location in the hierarchy. Similarly, in Newcomb's Fraternity, friendships with highly-ranked brothers may confer greater prestige than those with lower-ranked ones.
In addition to the likelihoods, we can also compare the memory estimateλ across models and data sets. Since the model assumes that the impact of past endorsements decays at rate λ, the quantity t 1/2 = − log(2)/ log(λ) represents the half-life of system information according to the inferred dynamics, in units of observation periods. In the Math PhD data, the favored Root-Degree score gave a half-life of t 1/2 ≈ 5 years. In the Parakeets data, the half-life estimated under SpringRank is t 1/2 ≈ 1.7 weeks for the first group and t 1/2 ≈ 0.8 weeks for the second. The small number of observation periods implies that these estimates should be approached with caution. Finally, in the Newcomb Fraternity data, the SpringRank half-life was t 1/2 ≈ 2 weeks. This suggests that the friendships in this data set evolved on timescales much shorter than the full semester. This likely reflects the fact that the brothers did not know each other prior to data collection, requiring them to form their social relationships from scratch. An important caveat in interpreting these estimated half-lives is that the indirect influence of an interaction may extend far beyond its direct influence. In the Math PhD data, for instance, while the half-life indicates that only a quarter of hiring events will be directly "remembered" in the system after a decade, those events will have influenced ten cycles of hiring, which may further reinforce the patterns established by the earlier events.

Root-Degree PageRank SpringRank
Math  Figure 1) compared to the mean critical value β c 1 for each system. β c 1 is calculated as in Theorem 1, using as m the mean number of interactions per time-step in the observed data. As in Table 1, the parameters corresponding to the highest log-likelihood are shown in bold. Estimates shown with an upper asterisk * exceed the approximate critical value by two standard errors, while estimates shown with a lower asterisk * are smaller than the approximate critical value by two standard errors. See Fig. 9 for simulated trajectories using the inferred parameters.
As described in Theorem 1, in the long-memory limit, our model has distinct egalitarian and hierarchical regimes, separated by a critical value β c 1 . The model's estimate of β 1 allows us to roughly locate empirical systems within these regimes. There are two necessary points of caution. First, when the estimateλ is far from the idealized long-memory limit, hierarchical and egalitarian regimes may not be sharply distinguished. Second, in the Math PhD and Parakeet data, the number of updates m varies between time steps. Here, a reasonable approximation is to use the average number of updatesm per time step. Using this average and Theorem 1, we computed an approximate long-memory critical value β c 1 for each empirical system.
Comparing the data-derived preference estimatesβ 1 to the approximate critical values β c 1 reveals that all four empirical systems are in or near the hierarchical regime ( Table 2). The Root-Degree estimates of β 1 tend to be very close to the approximate critical point. For the Math PhD data, in which Root-Degree is the preferred model, the estimate of β 1 is slightly but statistically-significantly below the critical value. In each of the other three data sets, the estimate is slightly above the critical value, and significantly so in the two Parakeet groups. Given the presence of a bistable regime in the Root-Degree model (Fig. 3(a)), the estimate of β 1 for the Math PhD data is consistent with persistent hierarchy despite the fact that the estimate falls slightly below the critical threshold. Indeed, simulations with the inferred parameters produce persistent hierarchical structure similar to that observed in the data (SI Appendix, Fig. S5). The PageRank estimates behave similarly to Root-Degree, although the finding in Parakeets (G2) is reversed. The presence of a bistable regime in the PageRank model ( Fig. 3(b)) indicates that these findings are consistent with persistent hierarchy in any of these data sets (see SI Appendix, Fig. S5 for simulated dynamics). Finally, in the SpringRank model, which obtains the highest likelihood for both Parakeet data sets and the Newcomb Fraternity, the estimated values of β 1 significantly exceed the estimated critical values, and tend to lie in or near the range [2,3]. In summary, all three models suggest that the system corresponding to each data set is in or near the regime of self-reinforcing hierarchy. Our model also assigns intepretable, time-dependent ranks to empirical data (Fig. 4). For the Math PhD Exchange network, for example, the raw placement share (Fig. 4a) and Root-Degree model (Fig. 4b) show strong qualitative agreement, with institutions that place the most candidates occupying higher ranks. Due to the relatively large estimatesλ, both the Root-Degree and PageRank models (Fig. 2b-c) produce smoother rank trajectories than the purely-descriptive placement share with 8-year rolling average. In contrast, the SpringRank score generates qualitatively different trajectories that are less sensitive to raw volume (Fig. 4d). For instance, SpringRank places Harvard at the top over most of the time period, while the other scores prefer MIT. This difference reflects SpringRank's sensitivity to where Harvard's graduates were placed, a consideration which Root-Degree entirely ignores. Similarly, SpringRank places Chicago and Yale noticeably higher than Wisconsin-Madison, despite all three having similar numbers of placements.

discussion
We have proposed a simple and flexible model of persistent hierarchy as an emergent feature of networked endorsements with feedback. When the preference for high status exceeds a critical value, egalitarian states destabilize and hierarchies emerge. The location of this transition depends on the structure of the score function and of the node's preferences. Our findings emphasize that winner effects do not require internal, rank-enhancing feedback mechanisms. Social reinforcement through prestige preference is sufficient to generate social hierarchies. Crucially, our model has a tractable likelihood function, supporting principled statistical inference of parameters-for both preferences and memory strengthfrom empirical data. In the four data sets analyzed, we found that links are typically formed in alignment with the hierarchy (β 1 > 0) but that they are preferentially created to other nodes with similar ranks (β 2 < 0). The likelihood also opens the door to model selection to determine relevant score functions. We found that networked ranking methods that capture transferable prestige are preferred over non-networked methods in some but not all systems. Due to its flexibility, our framework can be applied to additional data sets, score functions, and/or preference models to test the generality of these empirical observations.
There are limitations to our approach. First, we specified a fixed parametric form for the utilities with Eq. (4) and Gumbel-distributed noise with Eq. (5). Other choices may be more justified in particular applications, ideally informed by domain-specific considerations. Importantly, our inferential framework allows for quantitative evaluation and comparison of these choices. Taking advantage of this, future work could systematically explore the most appropriate functional forms in systems from diverse scientific domains. Second, our model assumes that all nodes use identical preference parameters β 1 , β 2 and score vector s when computing utilities. The latter is an especially strong assumption, since it requires each node to have global knowledge of the endorsement network, or at least of the score vector. This is unlikely to be true in real systems, and should be regarded as a modeling device. Future work, along the lines of Hobson and DeDeo [2015], could explore the interplay between the cognitive capabilities of individuals represented by nodes and the information available to them in the formation of social hierarchies.
Our model points to several other avenues for further work. A crucial step would be to extend extant network-based models [König and Tessone, 2011;König et al., 2014;Krause et al., 2013] so that their parameters could be statistically learned from data. This would enable comparative validation of different modeling frameworks. Studies of the relationship between measures of time-dependent centralities [Taylor et al., 2017[Taylor et al., , 2019Liao et al., 2017] and dynamic models of hierarchy would also be valuable. In particular, the theory of time-dependent centralities faces an important methodological issue: different reasonable ranking methods can yield directionally different orderings of nodes when applied to the same data set [Mariani and Lü, 2020]. Their performance on external validation tasks, such as the prediction of central nodes in spreading processes [Lü et al., 2016], may also vary significantly. Because the theories of centrality and generative networks have evolved largely separately, evaluating the suitability of a centrality metric for a given dynamic system can be difficult. Our inferential approach offers a candidate validation task to overcome this challenge: good centrality metrics are those which most effectively predict the future evolution of the system. This approach enables us to not only compare different score and utility functions in a principled manner but also explore their relative importance in observed networks. For instance, one could study the relative influence of degree-based and SpringRank scores by incorporating both into our model and then analyzing their distinct coefficients. Further work in this direction could reveal how different forms of centrality combine to govern the evolution of interaction networks. We anticipate that a fruitful dialogue between centrality theory and generative models of time-varying networks will deepen our understanding of the feedback mechanism between local interactions and hierarchical structures.

published article
This article was published in the Proceedings of the National Academy of Sciences. We request that all citations refer to the published version. software A repository containing our data, model implementation, and figure generation scripts is available at https://github.com/PhilChodrow/prestige _ reinforcement.

acknowledgments
We thank the anonymous reviewers for their constructive feedback. We also thank Dakota S. Murray (Indiana University, Bloomington), Kate Wootton (Swedish University of Agricultural Sciences; University of Colorado Boulder), VPS Ritwika (University of California, Merced), and Rodrigo Migueles Ramírez (McGill University) for extremely helpful discussions during the early phase of this work. We are also grateful to the organizers of the Complex Networks Winter Workshop in Quebec City at which this work was conceived, and to the inhabitants of Café Félin Ma Langue Aux Chats, Quebec City, Canada for contributing to an exceptionally productive working environment. MK acknowledges support from the Army Research Office Grant W911NF-18-1-0325. PSC acknowledges support from the National Science Foundation Award 1122374. DBL acknowledges support from the National Science Foundation Award SMA 1633791 and the Air Force Office of Scientific Research Award FA9550-19-1-0329.

gender representation in cited work
Recent work in several fields of science has identified gender bias in citation practices: papers by women and other minoritized groups are systematically under-cited in their fields [Mitchell et al., 2013;Dion et al., 2018;Caplar et al., 2017;Maliniak et al., 2013;Dworkin et al., 2020]. Gender bias can arise through explicit and implicit bias against a person's known gender identity, or against a name commonly used by members of a marginalized gender identity [MacNell et al., 2015;Paludi and Strayer, 1985;Moss-Racusin et al., 2012]. In this work, we proactively sought to include relevant citations by non-male authors.
We manually gender-coded the authors in the works cited according to personal acquaintance, instances of pronoun usage online, or first name. This method is limited: names and pronouns may not be indicative of gender; gender may change over time; and manual coding is inherently subjective and subject to error. We focused on the first and last authors because typically, though not always, the former is the leading researcher and the latter the senior author in the disciplines included in our references. Of the works cited in the main text, 29% had a non-male first author and 13% had a non-male last author. Of those with at least two authors, 38% had either a non-male first author or a non-male last author.
This statement is modeled on those found in Torres et al. [2020] and Dworkin et al. [2020]. We join those authors and many others in calling for collective effort to promote equitable practices in science.
references [Ball and Newman, 2013]   In this section, we prove a set of linear stability results that generalize Theorem 1 in the main text. Our generalizations account for (a) nonlinear features and (b) multiple updates per round.
Throughout this section, we consider a utility function of the form where each φ : R n → R n×n is a smooth feature map; β ∈ R is a preference parameter indicating relative importance of the th feature; and φ ij (s) is the ijth entry of φ (s). We collect the parameters β in a vector β ∈ R k . The utility function in Eq. (4) from the main text is a special case with linear feature map φ 1 ij (s) = s j , and quadratic feature map, φ 2 ij (s) = (s i − s j ) 2 . We also define the rate matrix G = [n −1 p ij ], whose (i, j)th entry gives the probability that, in a given time step, node i chosen uniformly at random endorses node j (see Eq. (5) in the main text for the definition of p ij ).
Since we aim to characterize the linear stability of egalitarian fixed points, we will consider the Jacobian of the rank vector γ evaluated at egalitarian fixed points. We will therefore evaluate the Jacobian at s 0 = θe, where θ ∈ R. By definition, γ = n −1 G T e = n −1 i γ i , where γ i is the ith column of G. Differentiating and applying the chain rule, we have where Γ i = diagγ i and φ i· (s 0 ) is the ith row of the th feature map evaluated at s 0 . At s 0 = θe, G = n −1 E. It follows that γ i = n −1 e and Γ i = n −1 I. We thus have We will express our primary results in terms of this matrix. When writing proofs involving dynamics, we will typically repress the timeargument of quantities like s and A. When time step t is implied, we will use the somewhat informal notation δs = s(t + 1) − s(t) and δA = A(t + 1) − A(t) to denote the increments of these and other quantities in the current time step.

A.1 DEGREE SCORES
Theorem 2 (Stable Egalitarianism with Degree Scores). When σ(A) = s = A T e, the vector s 0 = de is a root of f , where d = m n , and is the only egalitarian root. Furthermore, s 0 is linearly stable in the long-memory limit if and only if M(s 0 ; β) has eigenvalues strictly smaller than 1 m . Proof. We first derive the functional form of f . We can write Inserting this expression into (6), and recognizing n −1 G(t)e = γ(t), we have We can now check that s 0 is indeed the unique egalitarian root of f . Suppose that s = se for some scalar s. Then, which is only equal to zero when s = m n , as needed. Now computing derivatives, we have This matrix has strictly negative eigenvalues provided that the eigenvalues of M(s 0 ; β) are strictly smaller than 1 m , completing the proof. Corrolary 1. Using the Root-Degree score function, s 0 = m n e is a linearly stable fixed point of f if and only if β < 2 n m . Proof. It is convenient to treat the operation of taking the square root as part of the feature map, rather than part of the score function. We therefore suppose that s j is the in-degree of node j and that φ j (s) = √ s j . Computing from (8), we obtain M(s 0 ; β) = 1 2 This matrix again has a zero eigenvalue associated with the direction e. For any direction v ⊥ e, there is an eigenvalue 1 2 n −1 √ d β. From Theorem 2, s 0 will be linearly stable provided that as required.

A.2 PAGERANK SCORES
The PageRank score [Brin and Page, 1998;Page et al., 1999] is the solution s of the linear system where D o = diag(Ae). The Perron-Frobenius Theorem [Horn and Johnson, 2012] ensures that s is strictly positive entrywise. We assume s to be normalized so that s T e = n, which is contrary to the usual normalization s T e = 1. This choice amounts to a rescaling of the parameters β, and does not otherwise impact the analysis.
In the case of PageRank, it is difficult to derive a result for general features and we therefore work directly with the PageRank model with linear features. Proof. Uniqueness is a direct consequence of normalization: if s = θe and e T s = n, then we must have θ = 1. We next obtain a necessary condition describing roots of f . We start with a useful simplification. At any fixed point of f , we must have D o = mI. This is because, at any such fixed point, we must have A = mG, and nG is row-stochastic. For the purposes of analysis in the long-memory limit, we can therefore consider s to be defined by the simplified equation In the next time step, we will have αm −1 n(A T + δA T ) + (1 − α)n −1 E (s + δs) = s + δs .
Expanding and canceling yields The term o(1 − λ) includes terms involving the product (δA T )(δs), and relies on the fact that δs is a smooth function of A. Rearranging and dropping the asymptotic term, we obtain, in the long memory limit, This expression gives an implicit representation of f via the relation f (s, A) = lim λ→1

E[δs]
1−λ . We can therefore enforce f (s, A) = 0 by setting E[δs] = 0, obtaining the necessary condition E[δA T ]s = 0 for roots of f . Expanding this condition yields, Inserting (10) and rearranging yields the nonlinear system The largest eigenvalue of the matrix on the lefthand side is α −1 n −1 . This allows us to numerically solve (12) iteratively, by alternating between solving for s via a standard eigenvalue solver and updating G with the new value of s. This is the method implemented in the accompanying software and used to generate equilibria in Fig. 3. In order to derive the linear stability criterion, we divide both sides of (11) by 1 − λ and differentiate with respect to s, obtaining After inserting (10) and simplifying, we have The second line follows from the normalization of s, which implies that Es = ne, a constant vector which does not depend on s. Differentiating the righthand side then yields Evaluated at the egalitarian solution s 0 = e, this becomes To complete the argument, we note that, at the egalitarian solution of our model dynamics, A = n −2 E. Inserting and simplifying, we have Provided that α < 1, the premultiplying matrix on the lefthand side is invertible, and I − αm −1 n −1 E −1 = I + α(m − α) −1 n −1 E. This matrix has a single eigenvalue 1 + α(m − α) −1 with eigenvector e, and additional eigenvalues equal to unity in orthogonal directions. We then have In the PageRank-Linear model, M(s 0 ; β) = βn −1 (I − n −1 E), and we therefore have We can now read off the eigenvalues of J(s 0 ) analytically. The eigenvector e has eigenvalue −1, while any vector orthogonal to e has eigenvalue αβ − 1. This latter eigenvalue is strictly negative if and only if β < 1 α , as was to be shown.

A.3 SPRINGRANK SCORES
We return to the general formalism of score functions and features introduced at the beginning of this section. A SpringRank vector s for a matrix A with regularization α ∈ R is a solution to the linear system where, is invertible and s is therefore unique. Thus, throughout this section we will assume that α > 0, and correspondingly refer to s as "the" SpringRank vector of A. It is convenient to define L α = D i + D o − (A + A T ) + αI and Λ = D i − D o , in which case the SpringRank relation reads L α s = Λe.
Theorem 4 (Stable Egalitarianism with SpringRank Scores). When σ is the SpringRank map, the vector s 0 = 0 is a fixed point of f , and is the only egalitarian fixed point of the dynamics. This fixed point is linearly stable in the long-memory limit if and only if the matrix has eigenvalues strictly smaller than αn m . We will break the proof into a series of three lemmas. The first lemma calculates the analytical form of f . The second shows that s 0 = 0 is the unique egalitarian fixed point of the long-memory limiting dynamics f . The third gives the criterion for linear stability.
Lemma 1. The deterministic approximant f for the SpringRank vector is given by where L G = Γ + n −1 I − (G + G T ).
Proof. Let us fix an implicit time step t. Here and below, we use the notational template δM = M (t + 1) − M (t) to refer to increments in various quantities under the dynamics (1). For example, δA = A(t + 1) − A(t) refers to the increment in A under the dynamics. We compute directly We can also explicitly write out formulae for the increments in L α and Λ: where we have given a name to the Laplacian L ∆ = diag(∆ T e + ∆e) − ∆ T − ∆ of ∆. Note that δL α does not depend on α, and we therefore simply write δL = δL α .
We can now formulate a simple condition for equilibrium in expectation. We have Subtracting the SpringRank relation L α s = Λe from each side of this expression, we obtain (L α + δL)δs = (δΛ)e − (δL)s .
Since δL = O(1 − λ), the lefthand matrix is invertible in for small λ provided that α > 0. We therefore obtain . The term O((1 − λ) 2 ) arises from the product of O(1 − λ) and the copy of (λ − 1) within δΛ and δL. Taking expectations, . We next insert the expressions (15) and (16) and use the fact that E[∆] = mG. This gives We can simplify this expression by recalling that (L + αI)s = Λe by definition, as well as the identities Ge = n −1 e and G T e = γ. Inserting these identities and simplifying yields We now construct f , obtaining Since E[δs] = E[σ(λA + (1 − λ)∆)], we can write as was to be shown.
Lemma 2. When σ is the SpringRank map, the vector s 0 = 0 is a root of f , and is the only egalitarian fixed point.
Proof. To show that s 0 = 0 is a fixed point of f , it suffices to insert this solution into (14) and simplify, noting that, when s = 0, γ = n −1 e. To show that it is the unique egalitarian root realizable as a SpringRank score, suppose that se were a SpringRank score for some s = 0. Inserting this into (3) and using the fact that e is a zero eigenvector of the unregularized Laplacian, we would have The total in-degree must equal the total out-degree. Pre-multiplying by e therefore zeros out the righthand, leaving: which is a contradiction unless s = 0.
Lemma 3. The egalitarian root s = 0 is a linearly stable root of the SpringRank dynamics in the long-memory limit if and only if the matrix has eigenvalues strictly smaller than α m . Proof. We need to compute J(s 0 ), the Jacobian matrix of f at s 0 = 0. The fixed point will be stable provided that J(s 0 ) has strictly negative eigenvalues. To compute this Jacobian, we compute derivatives in (14). Doing so and applying the product rule, we have We calculate ∂L G ∂s in Equation (17) Evaluating this expression at s = 0, we have where L G must also be evaluated at s = 0. We have G(0) = n −1 E, which implies L G = 2(I − n −1 E). We insert this expression and the formula for ∂γ ∂s given in (8), obtaining Since L α is symmetric and positive-definite, L −1 α is as well. The stability of the egalitarian fixed point is therefore determined by the eigenvalues of the matrix inside the brackets. Multiplying by nm −1 , we find that a necessary and sufficient condition is that the matrix have eigenvalues no larger than α m , completing the proof.
Proof. We first compute the derivatives ∂(Gs) ∂s and ∂(G T s)

∂s
. The ith component of Gs is v i = j γ j s j . The product rule for scalar functions of vectors gives the ith row of the derivative: ∂Gs i ∂s = j γ j e j + j s j ∂γ j ∂s = γ + j s j ∂γ j ∂s .
Written in matrix notation, the first term is G. To write the second term in matrix form, note that we need to multiply ∂γ ∂s by the matrix each of whose columns is a copy of s. This matrix is S T . We therefore obtain ∂(Gs) ∂s = G + ∂γ ∂s S T .
To compute the second derivative, note that G T s = γ(e T s), with ith component γ i e T s. Using the product rule for scalar functions of vectors, we have ∂ ∂s γ i e T s = γ i e + (e T s) ∂γ i ∂s .
The first term will become the matrix whose ith row is γ i , i.e. G T . This yields ∂(G T s) ∂s = G T + (e T s) ∂γ ∂s .
Combining these expressions yields our formula for ∂L G s ∂s : as was to be shown.

b parameter estimation
Throughout this section, we use the shorthand {A(t)} = {A(t)} τ t=0 to refer to temporal sequences of matrices up to fixed time τ . We now describe a simple maximum-likelihood model for learning the parameter β from a sequence of observations {∆(t)}. By construction, ∆(τ ) depends on the sequence of state matrices {A(t)} only through the most recent state A(τ ). We may therefore factor the probability of observing the data given a set of undetermined parameters as: P({∆(t)}; A(0), λ, β) = τ t=0 P(∆(t); A(t), β) .
While the parameter λ has disappeared from the righthand side, this expression is nevertheless implicitly a function of λ since the value of A(τ ) given A(τ − 1) and ∆(τ − 1) depends on λ.
Let us write out a typical factor on the righthand side. Let k i = ∆ i· , and let K i = e T k i . Then, Taking logarithms and collecting terms that do not depend on β or λ into a constant C(t), we obtain log P(∆(t); A(t), β) = n i=1 n j=1 k ij (t) log γ ij (t) + C(t).
Standard theory of maximum likelihood in exponential families implies that L is convex in β for any fixed λ. This implies that, whenλ is known, we can solve forβ via standard first-or second-order optimization methods. Let L * (λ; {∆(t)}, A(0)) be the optimized loglikelihood for fixed λ. We then complete the maximum likelihood scheme by optimizing L * with respect to λ, which our accompanying software does via a customized hill-climbing algorithm. In general, L * may fail to be convex as a function of λ, and we therefore perform multiple runs with different initial values of λ in order to find the global maximum.
c additional model traces    Fig. 2). The parameters β 1 and β 2 are allowed to vary. Higher variances correspond to more strongly hierarchical states.