A scaling theory of armed conflict avalanches

Armed conflict data display scaling and universal dynamics in both social and physical properties like fatalities and geographic extent. We propose a randomly branching, armed-conflict model that relates multiple properties to one another in a way consistent with data. The model incorporates a fractal lattice on which conflict spreads, uniform dynamics driving conflict growth, and regional virulence that modulates local conflict intensity. The quantitative constraints on scaling and universal dynamics we use to develop our minimal model serve more generally as a set of constraints for other models for armed conflict dynamics. We show how this approach akin to thermodynamics imparts mechanistic intuition and unifies multiple conflict properties, giving insight into causation, prediction, and intervention timing.


I. INTRODUCTION
The battlefield is a scene of constant chaos.

-Napoléon Bonaparte
The unpredictability of armed conflict is cited in the classic texts on warfare, Sun-Tzu's The Art of War, Lanchester's Aircraft in Warfare, and Von Clausewitz's Vom Kriege. Despite seeming chaos in the midst of a single conflict, the ensemble of many conflicts displays multiple mathematical regularities including Richardson's law, the scale-free distribution of fatalities in interstate warfare [1,2]. How Richardson's law and other scaling patterns relate to one another remains unknown [3][4][5][6], but a framework unifying these and other conflict aspects could facilitate prediction or reveal hidden and spurious causes of surprising outcomes.
We show, by studying the Armed Conflict Location & Event Data (ACLED) Project [7], multiple quantitative regularities that we unify in a simple scaling framework [8]. Such regularities are evocative of scaling laws that emerge in disordered, driven physical systems [9], in animal societies with long temporal correlations in conflict dynamics [10], elections [11], cities [12], amongst other social systems [13]. We find that law-like behavior at sufficiently long scales in armed conflict data are captured by a randomly branching, armed conflict (RBAC) model. This model has an underlying fractal geography on which conflict "contagion" spreads, uniform dynamics of conflict development on this geography, and scale-free fluctuations in virulence, or intensity, between conflicts.
We extract these regularities from the statistics of conflict avalanches, consisting of spatiotemporally proximate events that have been joined into clusters. The clustered events consist of individual, localized conflict reports in ACLED, which serves as a database for conflict reports worldwide. Given that most of the data is from Africa, we focus on that region. Each conflict report is labeled by type of interaction, geographic location, date, estimated fatalities, involved actors, and other information (see Appendices of reference [8] for more details). Restricting our analysis to armed battles, we use a systematic definition for relating conflict events: we combine all conflict events within separation time a = 128 days and separation length b = 140 km to generate conflict avalanches. Thus, conflict avalanches define a set of spatiotemporally extended structures characterized by quantitative properties that, complementary to sociopolitical definitions of "battles" or "wars," are constructed only using physical measures of distance.
After having specified a and b (see reference [8] for further details), conflict avalanches can be described by total duration T , diameter L, infected geographic sites N (a measure of area), fatalities F , and number of conflict reports R. We discover that conflict properties display power law tails in distribution, scale nonlinearly with duration, and that the exponents characterizing both types of scaling are consistent with one another according to a minimal scaling hypothesis. Over the course of a single avalanche, each of these quantities increases monotonically with time. When they are averaged to generate the dynamical trajectories l(t), n(t), f (t), and r(t), we find that they are invariant under rescaling of the separation time. Taken together, these properties constitute phenomenological scaling variables describing how conflict starts from some epicenter, spreads in space and time, and generates conflict events like fatalities at infected conflict sites such as population centers. With this description as represented in Figure 1 in mind, conflict avalanches are reminiscent of cascades in other contexts like epidemiology [14], neural activity [15][16][17], or stress avalanches in materials [9], where universality and scaling provide valid, reduced descriptions of system dynamics. Despite tremendous social, cultural, and ecological complexity, the notion that conflict dynamics likewise conform to a similarly sparse description of conflict contagion is not only an intuitive analogy but one supported by quantitative evidence. We develop the model in Section II, building on pre- vious observations of how conflicts grow to motivate the model [18][19][20][21]. We show that the model is consistent with features of the data like functional forms, power law scaling, and exponent relations. For the reader's ease, we provide a table of all the variables discussed in this paper in Appendix Table S.1 and their estimated values from  data, model, and simulation are compiled in Tables I and  II. In Section III, we discuss the structure of the model and how it posits a minimal framework for conflict dynamics. Finally, we discuss insights for prediction and intervention in Section IV.

A. Model dynamics for conflict spread
We first draw a qualitative outline of our RBAC model. Imagine a big, compact region of length Λ that is susceptible to conflict. If conflict breaks out at a central site x i , it "infects" neighboring sites through transportation and social networks, growing outwards from the nucleation site x 0 to cover a set of sites x ≡ {x i }, a conflict avalanche of diameter at most Λ. At each newly infected region (e.g., township, county, province), conflict becomes endemic, generating instability, news reports, and fatalities. Far from the nucleating site, however, conflict potency will be lower as the relevance of distant conflict decays and the density of infrastructure supporting it shrinks (e.g., transportation networks [22]). Finally, conflict avalanches are characterized by spatiotemporal variation such that some regions or epochs show much more activity, a kind of spatiotemporally embedded virulence encoded in the intensity of nucleating events. As we see in Figure 2, that different regions show strongly varying levels of conflict is empirical fact. Deserts, mountains, and oceans show sparse conflict, if any, but such variation might also result from weak government [21,23], technological variation [24], or historical friction between ethnic groups [25]. 1 These elements of geography, endemicity, and virulence define the multiple, parallel processes in our model for armed conflict.
At the core of our model is a randomly branching tree representing the spread of conflict sites at which conflict events occur on the approximately two-dimensional surface of the earth. The tree has branches of average length B k at generation k, each of which give birth to an average of Q branches when they reach their branching points with resulting fractal dimension δ n = 1 + log(Q)/ log(B) as in Figure 3. The increasingly distant branching points of the tree mimic the way road networks become sparse far from highly interconnected cores [22]. At each time step, a randomly chosen branch is extended by unit length, reflecting the addition of a new conflict site on which conflict reports begin to accumulate. As a result, the time it takes for a site to reproduce-that is, seed another conflict site and further extend the conflict avalanche periphery-increases as the tree becomes FIG. 3. Random variation of "Nice Trees" grown for T = 8,000 time steps with branching number Q = 3 and varying branch length ratios B [26]. For battles, B = 6.6. There is one conflict site per unit length.
larger in a way reminiscent of how battle fronts spread [20]. These simple dynamics mean that conflict site number grows linearly with time having set n to share units with t in our model. Likewise, average avalanche diameter is determined solely by the fractal dimension after sufficient time, Eq 2 also defines the dynamical exponent ζ, which is equivalent to fractal dimension δ n under these minimal single-site growth dynamics. This minimal model capturing geographic spread cannot explain how conflict multiplies at each new infected location as is measured by reports and fatalities. In fact, the measured spatial dimensions for fatalities and reports apparently exceed the dimension in which they live, d F > d R > 2, because of conflict recurrence in fixed areas (Table I). This implies that in order to capture growth in social dimensions of armed conflict, we must account for a separate set of dynamics evolving on top of the geographic lattice. On each site x i that is infected on day t 0 (x i ), conflict becomes endemic and a cascade of conflict events begins. A cascade on site x i generates conflict reports as well as fatalities, the cumulative numbers of which we track as r xi (t) and f xi (t). A phenomenological scaling model for reports at site x i is The particularly poor exception is diameter growth l(t)/L. Data points are few for the longest avalanches, but we find long avalanches saturate the maximum diameter suddenly. Inspecting these avalanches in detail, we find they tend to hit hard boundaries like coastlines and national borders they cannot surpass. In the case of the Tunisian and Libyan revolutions, the aggregation of which is included in the shown average for the longest conflicts, the population is largely confined to the coastline. This suggests for conflict avalanches commensurate with geographic or political boundaries, it is essential to account for such boundaries delimiting their maximum extent.
with an analogous equation for f xi (t). Eq 3 accounts for site virulence v r (x i ) modulating local magnitude, growth scaling with exponent 1 − γ r shared across all conflict sites, peripheral suppression for sites that start later characterized by exponent θ r ≥ 0, and a finite rate at all times, = 1. When the growth exponent is at its maximum value γ r = 1, the new event rate decays quickly, and event count is solely determined by virulence, start time, and peripheral suppression.
By accumulating over the entire extent of the conflict avalanche, we obtain We expect that at large scales report growth scale with time, a scaling relation that defines the dynamical exponent δ r /ζ. In order to proceed with the calculation, we assume that random fluctuations in site virulence v r (x i ) are uncorrelated with the time at which a site appeared and use a mean-field approximation averaging over conflict avalanche extent, assumptions we verify with data later. Then, Eq 4 only depends on temporal and spatial scales, where we have denoted V r (x) ≡ v r (x i ) , the expected virulence over a single avalanche x, and the typical number of sites l(t) δn . With a single site added at every time step, the probability that any randomly chosen conflict site was first infected at time t 0 is uniform and for sufficiently large t. Normalizing Eq 7 by R x ≡ r x (t = T ) to remove dependence on conflict region virulence, we average over x to obtain the universal scaling function This presents our first exponent relation for growth in reports using the definition in Eq 5, A similar relation holds for fatalities f . Taken together, Eqs 1-9 describe predictions of functional forms and exponent relations that we verify with data.

B. Verifying dynamical model on data
Using the conflict avalanches that we construct with the data from ACLED as discussed in reference [8], we check whether or not predictions about universality and self-consistent exponent relations are supported by the data.
As one test of our predicted scaling form for normalized trajectories in Eq 5, we construct conflict avalanches after rescaling separation time a → 2a. Under such a change, conflict avalanches will increase in size and duration, although in a way that leaves the normalized functional form unchanged. We show in Figure 4 over an order of TABLE I. Dynamical exponents measured from Battles data, calculated analytically for RBAC model, and estimated from simulation (K = 10 5 samples). See Figure 7 for scaling in simulation.
magnitude of rescaling in a, the normalized scaling form to be well-preserved, confirming our predictions in Eqs 1, 2, and 8 that the dynamical trajectories do not change under temporal rescaling.
As another test of the dynamical hypothesis, we compare normalized trajectories of short and long conflict avalanches in Figure 4. We find that these trajectories largely collapse onto a universal profile-though national and geographic boundaries have an outsize effect on diameter growth for the largest avalanches. Importantly, these "finite-size" effects are not prominent after we include avalanche of all sizes, suggesting that scaling of this average is less sensitive to variation in boundary effects across avalanche scales. Taking note of this difference, we take normalized trajectories averaged over avalanches of all durations to measure dynamical scaling exponents δ n /ζ, δ r /ζ, and δ f /ζ shown in Table I. Given these trajectories, we can immediately check if the model exponent δ n /ζ = 1 is close to the data δ n /ζ = 1.06 ± 0.05, which confirms RBAC does indeed imitate the averaged geographic spread of real conflict avalanches across conflict regions and durations.
To check the predicted dynamical exponent expressions for δ r /ζ and δ f /ζ like in Eq 9, we must measure the site growth exponent γ r and peripheral suppression exponent θ r . First, we consider how to measure γ r . It can be measured directly from observed conflict trajectories for each site as given Eq 3. Taking its logarithm, we can fit for some constant and for some value of γ r such that We leave inside A the unknown combination of random virulence and exponent θ r as we discuss in further detail in Appendix Section A. After constructing conflict sites by taking Voronoi regions inside a conflict avalanche, we estimate γ r = 0.7 ± 0.2 and γ f = 0.6 ± 0.3 (we show the distributions of the exponents in Figure S.1), values we then use to calculate θ r .
To measure the decay exponent θ r , we compute how total activity at a site decays when it starts later in the conflict avalanche by combining decay profiles over all the sites over different avalanches. For a single site, the profiles are where we have defined the normalized time at which the site was infected g(x i ) ≡ t 0 (x i )/T and have assumed that the correction to first-order scaling going as 1/T is small. 2 Taking the average over sites x i within an avalanche and over conflict avalanches x (denoted by a bar), Eq 12 describes an averaged conflict event density by the relative time g that has passed, peaking at g = 0 and sharply suppressed at g = 1. This particular scaling collapse provides a prediction of how the density of events per site progresses during the course of the avalanche. Using our estimates for γ f and γ r , we use Eq 12 to fit the exponents θ f = 0.2 ± 0.3 and θ r = 0.4 ± 0.3 with 90% bootstrapped confidence intervals shown in Table I (see Appendix A for measurement details). Importantly, the resulting curves align qualitatively with our predictions as plotted in Figures 5 and S.2: the data show an increase in the conflict event rate at sites occurring near the beginning of the avalanche, with strong suppression at the end substantially different from when θ r = 0. With this confirmation, we combine our measured exponents to obtain 1 − γ r − θ r + δ n /ζ ≈ 0.9, which is remarkably Scaling form predicted in Eq 12 aligns qualitatively with the data given measured γr = 0.74. We show bounds on θr corresponding to 90% bootstrapped confidence intervals as θ − r and θ + r and RBAC simulation (orange). Since for each solution of θr there are corresponding fit parameters from Eq B1, the bounding lines for θ + r and θ − r indicate variability in the shape of the curve but not vertical displacement.
close to the measured value of δ r /ζ = 1.06. Similarly, 3, compared to the best fit estimate from Figure 4, δ f /ζ = 0.96. Though both of these exponent relations are satisfied within bootstrapped confidence intervals, there is substantial uncertainty in exponent values for conflict site dynamics γ r , γ f , θ r , and θ f such that the predicted relations are loosely bounded between δ r /ζ ∈ [0, 1.4] and δ f /ζ ∈ [0, 1.8]. That the best fit exponents conform closely to our predicted relations, indeed much closer than the uncertainty suggested by confidence intervals, demonstrates that our formulation aligns well with the dominant features of armed conflict growth. Thus, we find our mean-field formulation of conflict site growth in the RBAC model accurately captures site evolution, peripheral suppression, and tightly satisfies self-consistent exponent relations.

C. Conflict virulence and extinction
By definition, a conflict avalanche ends when the rate at which new reports ∂ t r xi (t) are generated falls below some threshold as is set by our separation time a. Then, conflict extinction is determined by when the most prolific site at time t falls below rate threshold C, Given t and looking over sites with starting times t 0 (x i ), the rate is dominated by the two peaks at the endpoints with starting times t 0 (x 0 ) and t 0 (x T ). As a result, the fastest rate is determined by the relative magnitudes of the exponents θ r and γ r . Since γ r > θ r , the rate at the core dominates, and the threshold is met when Vr and V f display power law tails whose measured exponents satisfy self-consistent equations derived from the scaling hypothesis (p > 0.8 compared to standard significance threshold p = 0.1) [2].
A universal constant threshold C would imply that V r ∼ T γr . More generally, we might expect that larger conflicts are more difficult to observe because of the "fog of war" or if resources for observation are limited such that smaller events do not register as easily [27]. Though our rate threshold is fixed by the separation time, a fluctuating observation threshold could be effectively represented by rate threshold C fluctuating with duration such as When ∆ r > 0, the threshold increases with conflict duration and thus size, implying that observers are unable to resolve the smaller events unfolding in the conflict. 3 In this more general case, the rate threshold condition in Eq 13 implies where the exponent γ r describes the decay of conflict event rate at any particular conflict site and exponent ∆ r describes how the ability to resolve individual conflict events fluctuates with virulence. Similarly, we can construct an argument for fatalities, which likewise leads to a dynamical scaling prediction for conflict virulence of fatalities. This scaling relationship between virulence and duration links local dynamics of conflict growth with conflict avalanche termination, a global property. To take this further, we ask what happens if the distribution of conflict virulence were distributed in a scale-free way, 3 On the other hand, ∆r < 0 presents the unlikely possibility that observations become more detailed with increasingly larger conflicts. Such an unrealistic outcome would suggest that this intuitive explanation is flawed, but we find reassuringly the sensible bound ∆r ≥ 0 to be satisfied.
Fluctuations in V r would thus induce scaling in conflict duration determined by predicted exponent relation, In order to verify this hypothesis, we calculate the virulence for every site in conflict avalanches using our estimates for γ r and θ r . We show the resulting distributions in Figure 6 for V r and V f , which both are statistically consistent with having power law tails. From the distributions, we determine β r = 3.0 ± 0.3 and β f = 2.5 ± 0.4. As has been previously noted [8], the distribution of duration P (T ) also displays a power law tail with α = 2.44±0.13. Then comparing virulence with duration T , we estimate the dynamical scaling exponent γ r + ∆ r = 0.66 ± 0.02. Interestingly, this measurement means that ∆ r = 0 is consistent with the data, and that the report rate threshold does not necessarily depend on the intensity of observed conflict. Taking this seriously, we remove an additional parameter by setting ∆ r = 0. This is in contrast to the same calculation for fatalities, ∆ f + γ f = 1.32 ± 0.05, which implies ∆ f > 0.3 given the bound γ f ≤ 1 (see Table I). Such a result suggests that conflict resolution for fatalities fluctuates, a conclusion that aligns with the difficulty of estimating fatalities accurately [7,27]. Reassuringly, these exponents satisfy the predicted scaling relation in Eq 18, and conflict avalanche extinction aligns with a universal threshold in a way consistent with our a universal separation time scale. Thus, we show the way that we relate virulence and duration, derived from assumptions about scaling and our definition of conflict termination, lead to self-consistent relations satisfied by the data.

D. Scaling framework
Beyond the scaling of virulence with final conflict duration, the way that the remaining scaling variablesdiameter L, extent N , fatalities F , and reports R-grow with duration also imply additional power law distributions, These are not assumptions but are mathematical consequences of unifying the conclusions in previous sections, and these power laws hold in the data as described at further length in reference [8]. The new exponents in Eq 19 are determined by relating site dynamics with total magnitude of conflict avalanche properties after accounting for virulence. Using fatalities as an example, we define the exponent combination d F /z, Thus, a positive exponent combination γ f + ∆ f means avalanches grow larger than uniform site dynamics on a branching tree would imply, the excess scaling captured in our model by conflict-site correlations induced by virulence. We calculate from the entries of Table I the contribution of such virulence. Fatalities show strong effects of virulence revealed by the difference 1.0 ≤ d F /z − δ f /ζ ≤ 2.3, consistent with γ f + ∆ f ≈ 1.3, and implying ∆ f > 0. Correspondingly with reports, we find that the exponent γ r = 0.74 accounts for the difference 0.6 ≤ d R /z − δ r /ζ ≤ 1.5 such that ∆ r = 0, consistent with a fixed conflict termination threshold as noted earlier. Virulence, however, seems to play little to no role in the geographic spread of conflict, 0.2 ≤ γ n +∆ n ≤ 0.8 and −0.1 ≤ γ l + ∆ l ≤ 0.4. This observation aligns with our model assumption that virulence is primarily a feature of the social dimensions of conflict but not of geographic spread. By connecting the dynamics of conflict growth with the distributions of conflict scaling variables, we unify within a single mathematical model all of these properties and confirm our hypothesis that social growth results from a combination of geographic spread and conflict virulence.

E. Simulation
As a final check, we simulate the RBAC model. We find close agreement with scaling patterns in the data as shown in Figure 7 and Tables I and II (see Appendix Section C for further details about the simulation).

III. A MINIMAL MODEL?
Our approach relies on scaling, self-consistency, and simple dynamical hypotheses to build a minimal model that unifies both social and geographic characteristics of armed conflict. Yet, there are sufficiently many components that one might ask if the model is overparameterized. We argue in this section that our model represents a dramatic simplification of the full space of possibilities encompassing 7 scaling variables (i.e., duration, diameter, extent, reports, fatalities, and two types of virulence) and their trajectories. In principle, each of the scaling variables constitutes an independent degree of freedom with infinitely more degrees of freedom for the shape of growth trajectories and their distributions. To specify the functional form of the joint probability distribution relating every such degree of freedom to one another without an informative prior is difficult given the sparse and noisy data available. Instead, we posit a form for the decomposition of the joint probability that is tractable and empirically verifiable starting with assumptions about scaling.
As an example, consider the growth of armed conflict in duration t, diameter l, and extent n. In the most general possible scenario, we have arbitrarily complicated functions relating each pair of variables. However, under our scaling hypothesis, we restrict ourselves to only consid-  Table I (dashed black  lines) and are similar to scaling in data (blue). Measured dynamical scaling functions are shown after having removed the nonzero intercept at t = 0 averaged over conflict avalanches with duration T ≥ 4 days. For n(t), we also require N > 1 and for f (t) that F > 2 fatalities. (right column) Distributions of scaling variables with exponents listed in Table II align closely. Distributions for both data and RBAC are shown above lower cutoffs and their scales matched such that the lower cutoffs coincide.
ering power law forms that correspond to three separate exponents, or degrees of freedom. Under self-consistency and the absence of any additional scaling, the third exponent must be determined in terms of the other two, leading to the relationship n ∼ t δn/ζ as follows from in Eqs 1 and 2. Adding onto this, we assume single-site growth dynamics, which imposes equality of fractal dimension and dynamical exponent, δ n = ζ. Hence, with the case of geographic growth, the combination of scaling, selfconsistent exponents, and minimal dynamics compresses an arbitrary number of degrees of freedom into a single degree of freedom captured by the scaling exponent δ n /ζ that we measure from data (blue triangle in Figure 8).
Bringing reports and fatalities into the fold as we show in the leftmost panel of Figure 8, our model can be represented as a graph of dynamical scaling variables. In particular, averaged reports growth r(t) is a function of geographic spread, given by δ n /ζ, uniform site dynamics specified by θ r and γ r , and mean virulence V r (x). Thus, each aforementioned component contributes an additional exponent to r as indicated by the four incoming arrows. By traversing this sparse graph and taking the exponent relation corresponding to each edge, it is possible to relate every dynamical scaling variable with any other, but note the absence of redundant edges: we have avoided specifying any more edges than necessary to connect all the scaling variables. This dynamical description of conflict growth reduces the open-ended problem of fitting conflict data to specification of a few exponents-to be precise one for the set t, l, n and two for reports r x (t) and two for f x (t)-whose relationships align quantitatively with the data. The mean virulence V r , however, is unusual as is indicated by its red text color in Figure 8. Unlike the other scaling variables in black, it is quenched and so does not change as conflict progresses. Instead, the criterion for conflict extinction relates it to the total duration, linking dynamics with fluctuations in conflict avalanche size. Thus, virulence plays a special role in our theory, driving the intensity of conflict site growth in a uniform way within the context of a single conflict avalanche but displaying scale-free fluctuations across many separate conflict avalanches. This aspect is represented in the midddle panel of Figure 8 as the power law distribution of virulence P (V r ). With this assumption, we can calculate distributions of all remaining variables using the dynamical scaling relations and obtaining vast simplification. For example, we can construct the distribution of fatality virulence P (V f ) by using the dynamical scal- and thus a power law form for the distribution P (V f ). Taken together, these componentsuniform growth dynamics, scale-free fluctuations in virulence, and avalanche extinction below some threshold rate-compose a set of mathematical relationships between measurable conflict properties that sparsely relate the many aspects of conflict. Beyond our model, these scaling relations serve as constraints delimiting the set of conflict models that, if specifying many further microscopic details and proposed mechanisms for conflict propagation, must still hew to the regularities that we find in the data.

IV. DISCUSSION
That the complex tangle of armed conflict reveals strong regularities at large scales is truly remarkable. As one notable example that might have led us to anticipate the opposite, consider the conflict avalanche spanning Tunisia and Libya [8]. This outbreak of civil wars, which was part of the Arab Spring, clearly adheres to the geometry of the coastline given the density of population there. In contrast with other conflicts, this war began with the end of dictatorship and devolved into infighting amongst multiple militias seeking control over land, natural resources, and government [28]. Furthermore, it is difficult to refute the argument that geography plays a defining role in this conflict avalanche's spread. Yet, in the face of many such particulars, the statistics that emerge from the ensemble display highly regular, emergent properties such self-consistent power law scaling and universal dynamics. Here, we exploit these regularities, using them to organize and unify social and physical properties of armed conflict in a scaling framework captured by our RBAC model.
Both qualitative understanding of conflict causes and observed regularities in the data motivate our starting assumption that multiple features of armed conflict abide by simple scaling laws [1,[3][4][5][6]. Although some of these features like the distribution of conflict sites might reflect a process external to conflict dynamics such as socioeconomic variability [21], it remains an open question of how such statistical patterns emerge in the first place. One set of hypotheses revolves around the idea that conflict is an example of self-organized criticality (SOC) [29]. Roughly speaking, one might imagine that slow growth of social tension contrasted with relatively abrupt conflict resolution leads to scale-free features [18]. This is a debated hypothesis, but we observe that SOC models such as forest fire models neither abide closely to our measured scaling laws nor account for the full set of conflict features [8,18,30]. At the least, SOC models must incorporate heterogeneity in space and time, which is, as we find, a defining feature of armed conflicts. Some physical analogs of these features like quenched disorder [9,31], dissipation [32], or repetition on sites [16,33] have been considered in canonical models for criticality in nonequilibrium phenomena-though armed conflict suggests variations on these themes that may apply to social phenomena. More generally, the features we measure and the relations we establish between them in the RBAC model present a set of quantitative constraints that can be brought to bear on other models for armed conflict dynamics.
One constraint of particular note for conflict models results from our hypothesis that spatial scaling in armed conflict arises from the underlying geography on which it evolves [34,35]. As a way of capturing the fractal nature of conflict site density, we assume that conflict sites form a randomly branching tree. In this scenario, conflict features are determined by transportation networks, population density, and other social factors [36]. In intriguing alignment, some data suggest that the number of intersections of a road is characterized by a power law with exponent 2.2 ≤ u ≤ 2.4 [22]. Though conflict zones may be traversed in many ways, the overall statistics might be dominated by few major pathways such as the ring road in Afghanistan [35]. If so and if we think of intersections as meeting places where conflict actors converge, intersection density could account for why conflict extent is distributed with exponent u = 2.2. Further support for the idea that transportation networks influence conflict comes from results showing fractal dimension of metropolitan road networks globally span the range 1.2 ≤ D ≤ 1.7 [37], findings that are in agreement with our exponent for armed conflict extent δ n = 1.6. When a complete map of African transportation networks becomes available, it will be possible to further specify the mechanistic role of infrastructure on conflict dynamics.
Our approach reveals that conflict is not simply a geographic growth process but involves lattice-site dynamics resulting from its social nature. In particular, the density of reports and fatalities surpasses the two-dimensional physical landscape in which they are embedded, showing that the temporal dynamics at each lattice site are relevant. At each conflict site, reports and fatalities grow independently of geographic spread and are only rescaled in magnitude by final conflict duration. This suggests that conflict spreads locally in a common way-perhaps from shared social network structure across different parts of Africa or universal conflict spreading dynamics [38]. This would suggest that universality in conflict manifests in both local structure as well as in the statistics across many conflicts that span larger scales [39]. Overall, we find armed conflict dynamics are a consequence of underlying geography, asymmetry in between the core and periphery, and conflict virulence, aspects that are expressed through the scaling exponents. Interestingly, our model reveals the presence of correlated fluctuations in conflict intensity, or conflict virulence, indicating spatiotemporal disorder separate from universal dynamics. Virulence specifically enhances fluctuations in social dimensions, reports and fatalities, in our model (though exponent differences suggest that some analog of virulence, e.g., population density, may matter for spatial extent, its effects are much weaker). Superlinear scaling of social phenomena with population number has been observed in the dynamics of cities and has been argued to promote innovation and growth [12], but social scaling might likewise facilitate the spread of conflict, disinformation [40], or disease [41]. This aligns with the possibility that virulence reflects local social properties such as weak governance (e.g., comparing South Africa with Eastern Somalia [42]) or, similarly in primate societies, weak conflict management by leaders [43,44]. Alternatively, virulence could reflect a property of the instigating set of events as in primate society in which conflict duration grows with originating event severity [10]. Importantly, our finding of correlations in intensity over time suggests final conflict properties might be predicted at the onset. That conflict extinction is determined by the rate of events at the core is consistent with scaling in the data, suggesting that the origin of conflict outbreak is quantitatively, and perhaps plainly, linked to conflict duration.
Besides highlighting the importance of granular, highresolution, and accurate social data to further the study of armed conflict [45], our work demonstrates the power of a thermodynamical approach to revealing and accounting for regularities in a complex and noisy social system [5]. If, as our minimal model suggests, geographic and social characteristics drive the evolution of conflict, then universality and scaling we observe may arise from the in-tersection of human social dynamics, regional properties like geography, and social structure. Here, we describe how we measure the conflict site growth exponents γ r and γ f and the virulence V r and V f . We measure γ r by using the functional forms for site growth as in Eq 10. To estimate the fitting parameters, we parameterize the logarithm of the scaling form to minimize the sum of two terms: one to fit the beginning of conflict avalanches and the other to fit the end. With reports as an example, We constrain the sum 1 − γ r ≥ 0. Then, we follow an analogous procedure for γ f . The resulting distributions are shown in Figure S.1. Given the long tail we find, we use the medians as estimates of the exponents instead of the means. Then, we take our best estimates of γ r and θ r , as described in the main text, to calculate the virulence per site at the end of the conflict avalanche, t = T . The averages of these measurements over all sites within a conflict avalanche returns the average V r , which we show in Figure 6.
Appendix B: Measuring θr and θ f To measure the peripheral suppression exponents θ r and θ f , we use the average profile defined in Eq 12. We parameterize the fit to include a coefficient determining units e A and a small "average" correction e B . The objective function for reports is the minimization problem where the averaged profile for reports depends on g implicitly through the relative time at which site x i started in the pertinent conflict avalanche. The form for B ensures that it remain positive (or zero) and the second, cost term ensures that it remain small as is assumed in the derivation of the scaling form. The weighting terms σ g are the standard deviation of our measurements used to obtain the averaged profile r xi (T )/T θr+γr−1 such that the fit is more tightly constrained by the more precisely estimated points. Finally, we discretize the relative time g ∈ [0, 1] to intervals spaced out by 1/9 as shown in Figures 5 and S.2. We solve Eq B1 using standard optimization techniques [46]. This procedure yields our initial estimates for the peripheral suppression exponents for the data.
For estimating the same exponents θ r and θ f from the RBAC simulation, however, there are two additional considerations that we take into account to solve the objective function defined in Eq B1. First, we are able to obtain long conflict avalanches and the singularity at t 0 = 0 becomes important to consider. Indeed, if we fit the profile with the first point at relative time g = 0.056, the emerging singularity at g = 0 can substantially distort the measured value of θ r . On some test examples, we find that the point at g = 0.056 jumps anomalously and forces the fit to match the remaining points poorly, an indication that our coarse-graining of g into intervals of 1/9 provides insufficient resolution to estimate θ r accurately when avalanches are much longer than typical ones in the data. However, it is the case that far from g = 0, nents γr and γ f estimated from regression to conflict site growth curves. Given this wide distribution, we take our best estimate of the exponent to be the median with confidence intervals given by the 5th and 95th percentiles as given in Table II. the singularity has much smaller effect and by simply excluding the point at g = 0.056, we recover accurately θ r = 0.5, the value cited in the main text. Though in principle similar bias is also an issue for g = 1, it does not skew our estimate of the exponents strongly and so we include it to replicate the procedure we use for the data as closely as possible. The second modification we make to the fitting procedure comes from the fact that σ g is no longer dominated by sampling noise and reflects the fact that fluctuations become larger near the singularity at g = 0. Since fluctuations in the model are a function θ r , the objective behaves deterministically with θ r , θ r is driven to large values, and the objective minimized by simply compressing the scaling function to vanishingly small values. For fitting the model, we replace σ g with e A such that the objective is rescaled by the typical value across the profile. We find that this allows us to get much more reasonable estimates for θ r and θ f while accounting for the typical scale of the average profile. Importantly, we find that these procedures lead to close fits of the averaged profile over the values of g that we consider. Putting these pieces together, we find close agreement between the exponents estimated from the model and data, providing a way of confirming the validity of our fitting procedures using the model.

Appendix C: RBAC simulation
We start by growing a randomly branching tree of fractal dimension δ n = 1.6 (calculated from taking the ratio of the separately measured exponents δ n /ζ and 1/ζ) emanating from a single seed site. Here, we consider Q = 3 and produce an initial set of three branches with an average extension factor B = 6.6. At each branching point, each set of children branches have random length B k (1 + η), where η is a random number chosen uniformly in the interval [−σ η , σ η ], σ η < 1 such that branches vary in length about the mean with fluctuations that grow proportionally with the mean. Given the lengths, the angle at which the branches split are chosen such that no branches will intersect with any other branches for a tree of arbitrary size. Examples of such random trees are shown in Figure 3.
On every newly added site, report and fatality dynamics are instigated such that the total number of events grow according to Eq 3. We set site dynamical exponents to their best fits: γ r = 0.74, θ r = 0.43, γ f = 0.56, θ f = 0.23, with V r sampled from power law distribution with exponent β r = 3 and lower bound of V r,0 = 1 to avoid very small conflict avalanches dominated by finitesize effects. At each conflict site, we treat the total cumulative number of events to be a continuous function of the discrete number of time steps t 0 (x i ) as would be the case in the limit of large avalanches. 4 This gives us the trajectories per site r xi (t) and conflict avalanche evolution r x (t) as well as the corresponding trajectories for fatalities, f xi (t) and f x (t).
Conflict avalanches are run til they reach the threshold rate of events determined by the scaling relation in Eq 14.
To simulate this, we take the random sample for virulence V r as mentioned above. Given a fixed, universal conflict rate threshold (e.g., C = 2 −7 , or one event per 128 days), the simulation ends when the mean event rate at the core crosses the threshold Thus, conflict extinction is determined by the combination of our fixed threshold for conflict rate, conflict avalanche virulence, and the universal rate with which it decays. The results are shown in Figure 7. fractal dimension from ?(t)/?(T ) ∝ (l/L) δ ? δ ? /ζ dynamical exponent for ?(t)/?(T ) ∝ (t/T ) δ ? /ζ ∆ ? part of virulence dynamical scaling exponent V ? ∼ T γ ? +∆ ? correction to singularity in site growth ζ dynamical exponent for diameter profile l(t)/L ∝ (t/T ) 1/ζ θ ?
peripheral suppression exponent Λ maximum cutoff length for conflict region ν diameter distribution exponent P (L) ∼ L −ν ξ correlation length τ fatalities distribution exponent P (F ) ∼ F −τ τ reports distribution exponent P (R) ∼ R −τ B conflict tree branch extension ratio C reports threshold for conflict extinction d ? /z dynamical exponent for ? ∝ T d ? /z f , F fatalities l, L diameter in km n, N number of conflict sites Q conflict tree branching number r, R number of reports rx i number of reports for site xi t0(xi) time of first event at conflict site xi t, T duration in days u extent distribution exponent P (N ) ∼ N −u v ? (xi) virulence for site xi V ? (x) virulence averaged over sites i in conflict avalanche x V ? virulence averaged over different conflict avalanches x conflict avalanche with sites {xi} xi conflict site i in conflict avalanche x z dynamical exponent for length L ∝ T 1/z