Urban skylines from Schelling model
UUrban skylines from Schelling model
F. Gargiulo , Y. Gandica , T. Carletti
1. Department of Mathematics and Namur Center for Complex Systems - naXys,University of Namur, rempart de la Vierge 8, B 5000 Namur, Belgium
We propose a metapopulation version of the Schelling model where two kinds of agents relocatethemselves, with unconstrained destination, if their local fitness is lower than a tolerance threshold.We show that, for small values of the latter, the population redistributes highly heterogeneouslyamong the available places. The system thus stabilizes on these heterogeneous skylines after a longquasi-stationary transient period, during which the population remains in a well mixed phase.Varying the tolerance passing from large to small values, we identify three possible global regimes:microscopic clusters with local coexistence of both kinds of agents, macroscopic clusters with localcoexistence (soft segregation), macroscopic clusters with local segregation but homogeneous densities(hard segregation). The model is studied numerically and complemented with an analytical studyin the limit of extremely large node capacity.
I. INTRODUCTION
Modern societies are often faced to segregation; dictated by race, religion, social status or incomes differences,it represents a major issue whose outcome can range from social unrest, to riots and possibly to civil wars. Theunderstanding of the rise of such phenomenon has thus attracted a lot of attention from economists, politicians andsociologists [1–4].In a couple of papers written in the late 60s, Thomas Schelling [5, 6] proposed a stylized model to describe theunset of segregation and pointed out a counterintuitive but widely observed result: a well integrated society canevolve into a rather segregated one even if at individual level nobody strictly prefers this final outcome. Even whenindividuals are quite tolerant to neighbors of their opposite kind, allowing them to relocate themselves to satisfytheir preferences - namely maximize their perceived local fitness/utility - will make segregation to emerge as a globalaggregated phenomenon not directly foreseen from the individual choices.Since the pioneering works of Schelling the model has attracted the attention of the community of physicists andmathematicians, interested in its simplicity and in the emergent behaviors recalling models such as Ising and Pottsones [9–15]. In its simplest form the model proposed by Schelling considers a population composed by two kinds ofagents sitting on the top of a regular lattice (or a 1 dimensional ring), each site being able to receive at most one agent.The latter being allowed to hop to a new empty lattice site once unhappy, that is once the fraction of agents of heropposite kind in her Moore neighborhood, is larger than a given tolerance threshold . The striking result by Schellingis that segregation will emerge for tolerances slightly larger than 1 /
3, well below the “more natural bound”1 / II. THE MODEL
In this paper we elaborate further in this direction by considering a metapopulation version of the Schelling modelwhere two kinds of agents, say Red and Blue, can move across N nodes, each of which can receive at most L agents (carrying capacity) at any given time. We assume there are in the model ρLN vacancies, i.e. empty spaces, and anequal number of Red and Blue agents. To make things simple we assumed a constant carrying capacity for each node, but of course one can improve the mode by consideringa different value of L i for each node i . a r X i v : . [ phy s i c s . s o c - ph ] M a y The considered model resembles to the one proposed in [16], the main difference being that in our case agents don’thave any information about the selected destination, therefore also moves that decrease or leave invariant the fitnessare allowed, we define such case weak liquid version of the Schelling model, the liquid case referring in the literatureto agents’ moves for which the fitness does not decrease [9]. A second difference is that in [16] the agent fitness iscomputed using single node informations, namely the number of agents in a given sites, on the contrary we reintroduceas Schelling originally did, the concept of spatial proximity [6], agent fitness takes into account the number of agentsin a given node and in the neighboring ones, namely nodes at distance 1 from the current node (local fitness).Nodes are assumed to be arranged in regular lattices, as initially assumed by Schelling, but the model can be easilyextended to complex networks. The local update rules is defined as follows. The neighborhood of an agent is given bythe topological neighborhood of the node where she lives, including the latter. At each time step an agent is selectedand her fitness is computed as the fraction of agents of the opposite kind of her, living in her neighborhood withrespect to the total number of agents living in the same neighborhood. Mathematically, assuming she is a Blue agentliving in node i , then her fitness is given by: f Bi = (cid:80) j ∈ i n Aj (cid:80) j ∈ i ( n Bj + n Aj ) , (1)where n Xj , X = A, B , is the number of agent of X -kind in node j , and we used the notation j ∈ i to denote all nodes j belonging to a neighborhood of node i , including the latter, that is the set of nodes at distance smaller or equal to1 from i .As in the original Schelling model, agents are unhappy if their fitness is larger than a given tolerance threshold , (cid:15) ∈ (0 , k not completely full, i.e. n Ak + n Bk < L :if f Bi > (cid:15) ⇒ agent A leaves node i . (2)The case for a Red agent is similar.One time step is the random selection with reinsertion of (1 − ρ ) N L agents. We define the convergence time , to bethe time needed for the system to reach the equilibrium, namely once no agents will move anymore.
III. RESULTS
We hereby present the numerical analysis of the proposed model once the underlying network is a regular latticewith periodic boundary conditions and each node has 4 neighboring nodes. The system is initialized with ρN L vacancies, (1 − ρ ) N L/ N = 400 nodes. The carrying capacity has been fixed to L = 100 and we check that the initial conditions satisfy thelocal constraint n Ai + n Bi ≤ L for all i . Throughout the paper the emptiness has been fixed to ρ = 0 . A. Single node properties
The aim of this section is to present the local properties of the system, namely at the level of single nodes. Themetric we used is the average value , over all the nodes, of the node magnetization : (cid:104) µ (cid:105) = 1 N (cid:88) i | n Bi − n Ai | n Bi + n Ai , (3)small values of (cid:104) µ (cid:105) mean that, on average, each node is populated by the same number of agents of both kinds, whilelarge values are associated to nodes filled with agents of only one kind.For (cid:15) ≤ . (cid:104) µ (cid:105) →
1, meaning that thesystem stabilizes into a frozen state where local segregation is present in all nodes: each node contains only agentsof one kind. Let us observe (see Fig. 1 B upper plot) that the same behavior is also present in the simplified modelwhere the fitness is calculated on the single node, as done in [16], and thus it is intrinsic to the displacement dynamicsand not to the way the agent fitness is computed. We also notice (see Fig. 1 panels A and Fig. 6 that as (cid:15) decreasestoward zero, the time needed to reach the frozen state gets longer, going to infinity in the limit (cid:15) → quasi-stationary non-segregated state (see red circles and orange stars curves in Fig. 1A). time . . . . . h µ i ✏ =07 ✏ =06 ✏ =05 ✏ =04 ✏ =03 ✏ =028 . . . . . < µ > < n m a x > . . . . . . . . . ✏ . . . . . < f e m p t y > P ( n ) n " = 0 . " = 0 . " = 0 . " = 0 . " = 0 . A B C
SEGREGATED CELLS MIXED CELLS
FIG. 1: Single node behavior. Panel A: average magnetization (cid:104) µ (cid:105) = N (cid:80) i | n iA − n iB | n iA + n iB as a function of time for a single genericsimulations. For (cid:15) ≤ . i ( n Ai + n Bi ) (middle plot) and fraction of emptynodes f empty (lower plot) as a function of (cid:15) . For (cid:15) ≤ . (cid:15) → n = n A + n B ) for few values of (cid:15) = 0 . , . , . , . , .
7. Averages shown inPanels B and C have been obtained performing 100 replicas of the model.
A second fundamental self-organized phenomenon emerges for (cid:15) ≤ .
5, the initial homogeneously distributed popu-lation organizes itself into an heterogeneous state across the network nodes (Fig. 1 panel C lower plots): most of thenodes contain ∼
10 agents, while very few nodes have as much as ∼
100 agents, recall that L = 100 is the maximumnode capacity.Observe that such asymptotic distribution is correlated with the time the system spends in the quasi-stationarynon-segregated state, the longer this time the more the population distribution across nodes moves from a Poissoniandistribution (Fig.1 panel C upper plots) to a power law (Fig. 1 panel C lower plots). The maximal node population n max = max i ( n Ai + n Bi ) increases for (cid:15) → L = 100) is never reached. At the same time we observe the formation of a relevant fraction of completely emptynodes, i.e. nodes for which n Ai + n Bi = 0 (Fig. 1 panel B lower plot).The transition to the local segregation at (cid:15) = 0 . ∀ i n iA n iA + n iB ≤ (cid:15) n iB n iA + n iB ≤ (cid:15) , (4)but it is straightforward to observe that this equation admits a solution, for which n iA > n iB >
0, only for (cid:15) > /
2. Hence for (cid:15) ≤ / B. Global properties
To analyze the spatial structures, thus beyond the single node, we define a cluster based on node’s majority, moreprecisely two linked nodes belong to the same cluster if they both are characterized by the majority of agents of thesame kind : i, j ∈ C if { n Ai > n Bi and n Aj > n Bj } or { n Ai < n Bi and n Aj < n Bj } . (5)If a node share the same (positive) number of Red and Blue agents, then it will be considered part of the interface .An edge between two neighboring nodes i and j is considered interface if ( n Ai − n Bi )( n Aj − n Bj ) ≤
0. Hence a clusteris made by nodes while an interface can contains both nodes and edges among them.This indicator allows us to show that (Fig. 2 panels B and C) for (cid:15) ≤ . . N and 0 . N and at the same time theaverage sizes of the first and second largest cluster added together cover more than 75% of the available nodes. Thiscritical threshold is the same observed in the original Schelling model. Notice that for (cid:15) = 0 . (cid:104) S max (cid:105) for0 . < (cid:15) ≤ . soft segregation because itallows a small mixing in the population. For (cid:15) ≤ . hard segregation of the population. For (cid:15) → . < (cid:15) ≤ . t − /z ( z = 4 for (cid:15) = 0 . z = 3 for (cid:15) ≥ . (cid:15) ≤ . IV. AN ANALYTICAL SIMPLIFIED MODEL
To gain insight into the behavior of the previously presented model, we hereby introduce a model ables to capturethe main behavior of the metapopulation Schelling model, but simple enough to be analytically tractable. Because weconsider the case for extremely large L , our model complements the one proposed in [16] devoted to the case L = 2.Using the notations introduced before the state of the system is thus completely characterized by the knowledge of( (cid:126)n A ( t ) , (cid:126)n B ( t )), being (cid:126)n A ( t ) = ( n A ( t ) , . . . , n AN ( t )) and (cid:126)n B ( t ) = ( n B ( t ) , . . . , n BN ( t )). For a sake of simplicity we decidedto compute the fitness using only the information from a single node, Eq. (1) with j = i . The system evolution isdone as previously and we still use the weak liquid version, an unhappy agent will move to an uniformly randomlychosen new node, provide there is enough space there.The model is thus intrinsically stochastic and hence it can be described by the probability P ( (cid:126)n A , (cid:126)n B , t ) to be attime t in state ( (cid:126)n A , (cid:126)n B ), whose evolution is dictated by the master equation: P ( (cid:126)n A , (cid:126)n B , t + 1) = P ( (cid:126)n A , (cid:126)n B , t ) + (6)+ (cid:88) ( (cid:126)n A (cid:48) ,(cid:126)n B (cid:48) ) (cid:2) T ( (cid:126)n A , (cid:126)n B | (cid:126)n A (cid:48) , (cid:126)n B (cid:48) ) P ( (cid:126)n A (cid:48) , (cid:126)n B (cid:48) , t ) − T ( (cid:126)n A (cid:48) , (cid:126)n B (cid:48) | (cid:126)n A , (cid:126)n B ) P ( (cid:126)n A , (cid:126)n B , t ) (cid:3) , Excluding the case where there is not strict majority in a node, this definition can be restated as: i and j belong to the same cluster if( n Ai − n Bi )( n Aj − n Bj ) > time ⇠ ( t ) ✏ =0.7 ✏ =0.6 ✏ =0.5 ✏ =0.4 ✏ =0.3 . . . . . . . ✏ . . . . . . . < S , S > first clustersecond cluster s ✏ . . . . . P ( s ) . . . . . . ✏ = 0 . ✏ = 0 . ✏ = 0 . ✏ = 0 . ✏ = 0 . n A > n B n B > n A n B = n A = 0 n A = 1 , n B = 0 n A = 0 , n B = 1 n A = 0 , n B = 0 A1A2 A5A4A3B CDNO GHETTOSGHETTOS WITH COMPLETE SEGREGATION GHETTOS WITH PARTIALMIXURE
FIG. 2: Global behavior. Panels A: Nodes occupancy at equilibrium, each node has been colored according to the majority ofagents in it: white empty nodes, blue only B agents, R only red agents and shades for intermediate cases. Panel B: Cluster sizedistributions for several values of (cid:15) = 0 . , . , . , . , .
7. Panel C: First and second largest cluster as a function of (cid:15) . PanelD: Interface size as a function of time for several values of (cid:15) = 0 . , . , . , . , .
7. Results presented in panels B, C and D arebased on 100 replicas of the simulations. being T ( (cid:126)n A (cid:48) , (cid:126)n B (cid:48) | (cid:126)n A , (cid:126)n B ) the Transition probability to pass from state ( (cid:126)n A , (cid:126)n B ) to the new compatible one ( (cid:126)n A (cid:48) , (cid:126)n B (cid:48) ).For incompatible states we set T ( (cid:126)n A (cid:48) , (cid:126)n B (cid:48) | (cid:126)n A , (cid:126)n B ) = 0.The non zero transition probabilities can be computed using the rules previously defined (see AppendixB for adetailed account of the transition probabilities calculation), for instance the transition probability that an A agentmoves from node i –th to node j –th is given by: T ( n Ai − , n Aj + 1 | n Ai , n Aj ) == 1 N n Ai n Ai + n Bi Θ (cid:18) n Bi n Ai + n Bi − (cid:15) (cid:19) L − n Aj − n Bj L .
The function Θ( x ), defined to be 1 if x > X –agents in every node i at time t , (cid:104) n Xi (cid:105) ( t ) = (cid:80) (cid:126)n X n Xi P ( (cid:126)n A , (cid:126)n B , t ), X = A, B .Using the expression for the transition probabilities Eq. [B3] and assuming correlations can be neglected, i.e. (cid:104) ( n Ai ) (cid:105) ∼ (cid:104) n Ai (cid:105) , we obtain a system of finite differences describing the evolution of (cid:104) n Ai (cid:105) and (cid:104) n Bi (cid:105) .To go one step further, we introduce the average fraction of A and B in each node, α i = (cid:104) n Ai (cid:105) /L and β i = (cid:104) n Bi (cid:105) /L ,we rescale time by defining s = t/L and finally we assume each node to have an infinite large carrying capacity (seeAppendixB for a more detailed discussion): dα i ( s ) ds = − ρ α i α i + β i Θ (cid:18) β i α i + β i − (cid:15) (cid:19) + 1 − α i − β i N (cid:88) j α j α j + β j Θ (cid:18) β j α j + β j − (cid:15) (cid:19) (7)and dβ i ( s ) ds = − ρ β i α i + β i Θ (cid:18) α i α i + β i − (cid:15) (cid:19) + 1 − α i − β i N (cid:88) j β j α j + β j Θ (cid:18) α j α j + β j − (cid:15) (cid:19) (8)where ρ is the fraction empty nodes, (cid:80) i (1 − α i ( s ) − β i ( s )) = ρN .To disentangle the evolution of α i and β i we introduce a new set of variables, the node emptiness , γ i = 1 − α i − β i ,and the difference of fractions of A and B , ζ i = α i − β i . The new variables range in ζ i ∈ [ − ,
1] and γ i ∈ [0 , F ( x ) = x Θ(1 − x − (cid:15) ) + (1 − x )Θ( x − (cid:15) ) G ( x ) = x Θ(1 − x − (cid:15) ) − (1 − x )Θ( x − (cid:15) ) , we can rewrite Eqs. (B7) and (B8) as follows dγ i ( s ) ds = ρF (cid:18)
12 + ζ i − γ i ) (cid:19) − γ i N (cid:88) j F (cid:18)
12 + ζ j − γ j ) (cid:19) (9)and dζ i ( s ) ds = − ρG (cid:18)
12 + ζ i − γ i ) (cid:19) + γ i N (cid:88) j G (cid:18)
12 + ζ j − γ j ) (cid:19) . (10)The agents initialization used in the previous section, namely ρN density of vacancies and an equal density (1 − ρ ) N/ A and B agents, translate into ( ζ i , γ i ) uniformly distributed in a neighborhood of (0 , ρ ). We thus divide the domainof definition of ( ζ i , γ i ) into four zones (see Fig. 3) and we will look closely to the dynamics in the Z zone if (cid:15) < / Z zone if (cid:15) > / ! " -1+2 $ -1 1-2 $ Z Z Z $ ! " $ -1 -1+2 $ Z Z Z $ FIG. 3: The four zones Z i used to study the solutions of system (9), (10). Let (cid:15) < /
2, then one can easily prove that if ( ζ, γ ) ∈ Z one has F (cid:16) + ζ − γ ) (cid:17) = 1 and G (cid:16) + ζ − γ ) (cid:17) = ζ − γ ,so assuming that for all i one has ( ζ i (0) , γ i (0)) ∈ Z , then Eqs. (9) and (10) rewrite: (cid:40) dγ i ( s ) ds = ρ − γ idζ i ( s ) ds = − ρ ζ i − γ i + γ i N (cid:80) j ζ j − γ j , (11)as long as ( ζ i ( t ) , γ i ( t )) will not leave Z . Assume for a while this statement to hold, then the first equation can bestraightforwardly solved to give γ i ( t ) = ρ + e − t ( γ i (0) − ρ ), that is for all i , γ i ( t ) → ρ when t → ∞ . The secondequation can also be solved (see Eq. [S13]) and thus to prove that ζ i ( t ) → i when t →
0. This proves also aposteriori that ( ζ i ( t ) , γ i ( t )) will never leave Z .The average magnetization rewrites in such variables as: (cid:104) µ (cid:105) = 1 N (cid:88) i | ζ i | − γ i , (12)we have hence proved that for initial conditions in Z the magnetization asymptotically vanishes (see Fig. 4).The remaining case, (cid:15) > /
2, can be handle as well but it is more cumbersome (see AppendixC). Let us onlyobserve there that for ( ζ, γ ) ∈ Z one has F (cid:16) + ζ − γ ) (cid:17) = G (cid:16) + ζ − γ ) (cid:17) = 0, so assuming that for all i one has( ζ i (0) , γ i (0)) ∈ Z , then the right hand side of Eqs. (9) and (10) identically vanishes and thus the average magnetizationdepends on the domain where initial conditions have been set. However for a fixed size of the latter, one cannot satisfythe hypothesis ( ζ i (0) , γ i (0)) ∈ Z if (cid:15) is closer enough but larger then 0 .
5, indeed the Z zone shrinks to zero in thiscase (see Fig. 3). In this case one should take into account the dynamics of orbits whose initial conditions are set in Z and Z (see AppendixB and Figs. 8 and 9). In conclusion the analytical model perfectly fits with the ABM for (cid:15) ≥ . (cid:15) < .
5, because in this range the ABM dynamics is strongly dictatedby the stochasticity of the model, orbits tending to converge toward the equilibrium (0 , ρ ) ( (cid:104) µ (cid:105) ∼
0) are destabilizedby fluctuations and thus sent into the zones Z and Z ( (cid:104) µ (cid:105) ∼ . . . . . . . . . (cid:15) . . . . . . < µ > ABMmeanfield
FIG. 4: The average asymptotic magnetization as a function of (cid:15) . Each point is the average over 50 simulations whose initialconditions are close to the equilibrium point γ i = ρ and ζ i = 0. Black circles represent the results of the analytical model whilegrey diamonds to the ABM with 1–node fitness. Parameters are: N = 100 and ρ = 0 . V. CONCLUSIONS
We proposed and analyzed an extension of the classical Schelling model to a metapopulation framework, whose mainoutcome is the spontaneous emergence, for low values of the tolerance threshold, of heterogeneously populated nodeswithout any exogenous preferential attachment mechanism (for 1D lattice this phenomenon induces the formation ofurban skylines, Fig. 5). This behavior is connected to the permanence for long time of the system in a quasi-stationarynon segregated state, where in each node the two populations are equally distributed. This quasi-stationary statecan be recovered as stable equilibrium of the simplified analytical model. We can evince that the system stabilizationtoward the magnetized state (for the ABM) passes through the creation (by random agents moves) of highly populatednodes, hereby named towers . At the same time, global patterns emerge as in the classical Shelling model. Figure 5summarizes the possible behaviors of the system. For (cid:15) < . (cid:15) ≥ . (cid:15) . On the contraryfor low tolerance cases the clusters are formed around the towers that become stable points for a certain type of nodes(once an agent, whose kind corresponds to the majority already inside the tower, enters she never gets out). Oncea higher density zone starts to exist, this mechanism reinforces the (majority color) population growth in this nodeand in the neighborhood. The global patterns start therefore to stabilize around the towers.Local magnetization phenomenon has been explained using an analytical approach able to describe the differentequilibria of the model and the origin of the quasi-stationary state for low tolerances. NO GHETTOS NO LOCAL SEGREGATIONNO SKYLINES
SPATIAL GHETTOS NO LOCAL SEGREGATIONNO SKYLINES
SPATIAL GHETTOS LOCAL SEGREGATIONSKYLINES
SPATIAL GHETTOS LOCAL SEGREGATIONNO SKYLINES " = 0 . " = 0 . " = 0 . " = 0 . FIG. 5: Summary of the global outcomes of the Schelling metapopulation model on a 1D lattice. Each box represents a nodeof the 1-dimensional lattice; the height of the box is given by the number of Blue and Red agent inside the node. The coloredrectangle below the nodes represents the composition of the node obtained using the same color code used in Fig. 2 PanelsA1–A5.
Acknowledgments
The work of F.G., Y.G. and T.C. presents research results of the Belgian Network DYSCO (Dynamical Systems,Control, and Optimization), funded by the Interuniversity Attraction Poles Programme, initiated by the BelgianState, Science Policy Office.T.C. is grateful to Cyril Vargas, student at the ENSEEIHT - INP Toulouse (France), involved in a preliminary studyof this subject during his Master Degree. [1] Cutler D, Glaeser E (1997) Are Ghettos Good or Bad?. Quart. J. Econ. 112: 827-872.[2] Fagiolo G, Valente M, Vriend N (2007) Segregation in networks. J. Econ. Behav & org 64:316-336.[3] Frey W, Farley R (1996) Latino, asian and black segregation in U.S. Metropolitan areas: Are multiethnic metros different?Demography 33:35-50.[4] Conradt L, Krause J, Couzin I.D, Roper T.J (2009) Leading According to Need in Self-Organizing Groups. The Americannaturalist 173:305-312.[5] Schelling, T.C (1969) Models of Segregation. American Econ. Rev. 59: 488-493.[6] Schelling, T.C (1971) Dynamic models of segregation. J. Math. Soc. 1: 143-186.[7] Pancs R, Vriend N (2007) Schelling’s spatial proximity model of segregation revisited. J. of Public Economic 91:1-24.[8] Hatna E, Benenson I (2012) The Schelling Model of Ethnic Residential Dynamics: Beyond the Integrated - SegregatedDichotomy of Patterns. JASSS 15: 1-6.[9] Vinkovi´c D, Kirman A (2006) A physical analogue of the Schelling model. Proc Natl Acad Sci USA 103,51:19261-19265.[10] Dall’Asta L, Castellano C, Marsili M (2008) Statistical physics of the Schelling model of segregation. J. Stat. Mech7(L07002):1-10.[11] Schulze C (2005) Potts-like model for ghetto formation in multi-cultural societies. Int. Journal of Modernd Physics C16:351.[12] Stauffer D, Solomon S (2007) Ising, Schelling and self-organising segregation. Eur. Phys. J. B 57:473-479.[13] Grauwin S, Bertin E, Lemoy R, Jensen P (2009) Competition between collective and individual dynamics. Proc Natl AcadSci USA 106(49):20622-20626.[14] Rogers T, McKane A (2012) Jamming and pattern formation in models of segregation. Phys. Rev. E 85:041136 (1-5).[15] Rogers T, McKane A (2011) A unified framework for Schellings model of segregation, J. Stat. Mech 7(P07006): 1-17.[16] Durrett R, Yuan Z. (2014) Exact solution for a metapopulation version of Schellings model. Proc Natl Acad Sci USA111(39):14036-14041.
Appendix A: Patterns dynamics
The aim of this section is to present some details concerning the convergence of the agent based model to theasymptotic equilibrium pattern. We start by considering the convergence time , defined as the time needed for thesystem to reach such equilibrium where no agents has incentive to move anymore, and its dependence of the tolerancethreshold (cid:15) . Intuitively, if (cid:15) is large then most agents are happy and thus they will not move, hence the equilibrium willbe reached quite soon. On the other hand if (cid:15) is small, agents will be very often unhappy and thus they will relocatethemselves to reduce this uneasiness, this will increase the convergence time. In the limit (cid:15) → .
25 0 .
30 0 .
35 0 .
40 0 .
45 0 .
50 0 .
55 0 .
60 0 .
65 0 . ✏ = 0 . ✏ = 0 . ✏ = 0 . ✏ = 0 . ✏ = 0 . ⌧ "⌧ FIG. 6: Convergence time as a function of (cid:15) . The blue line represents the average convergence time as a function of (cid:15) for 100replicas of the system with same initial conditions and parameters values. The inset shows the boxplots for the convergencetimes for some representative values of (cid:15) = 0 . , . , . , . , . In Fig. 6 we report the results obtained by the simulation of the metapopulation model and we can observe that theconvergence time does not have a monotonic behavior as a function of (cid:15) . Passing from from (cid:15) = 1 to (cid:15) ∼ . (cid:15) correspond tothe absence of macroscopic clusters (see Fig. 5 main text). Let us notice that the time to converge increase between (cid:15) ∼ . (cid:15) ∼ . . ≤ (cid:15) < . (cid:15) = 0 .
5, this is thephase where macroscopic clusters are formed through a coarsening process (see Fig. 5 main text). Finally for (cid:15) ≤ . (cid:15) (see inset of Fig. 6).The patterns arising for small tolerance threshold are intriguing and peculiar to our model. As already observedthe system spends quite a long time in a quasi-stationary (almost) homogeneous state and then suddenly jumps toa macroscopic segregated one. Such behavior is schematically represented in Fig. 7 in the case of the 1 D -latticefor (cid:15) = 0 .
25. The system stabilization toward the magnetized state passes through the creation (by random agentsmoves) of a highly populated nodes, hereby named towers corresponding to monochromatic clusters in the 2D modelpresented before. Such towers become attracting selective places for each kind of agents: once an agent, whose kindcorresponds to the one of the majority already inside the tower, enters she never gets out from the tower because shewill be happy there and this will increase the unhappiness of agents of the opposite kind still inside the tower. Thenet result is that once a highly dense zone starts to exist, this mechanism reinforces the (majority color) populationgrowth in this node and in the neighborhood, because of the way the fitness is computed, using nodes at distance 1.0The global patterns start therefore to stabilize around the towers. As the first kind of agents stabilizes, with a certaindelay also the stabilization of the second population is reached. Once the towers/clusters are formed, the completeequilibrium is reached as the interface between Blue and Red zones becomes empty. The lower is the tolerance value (cid:15) , the higher is the initial density unbalance needed to start the stabilization process. t=25000 t=20000 t=0 time . . . . . . . . . h µ i t=30000 t=35000 t=39500 . . . . . . . . . S , S FIG. 7: Generic time evolution toward a hard segregated pattern. Upper main plot: first and second cluster size as a functionof time. Lower main plot: average magnetization of the system as a function of time. Left and right panels: System snapshotsat different times. Each box represents a node of the 1-dimensional lattice, the height of the box is given by the total blue andred population of the node. The colored rectangle under the cell represents the majority color of the cell.
Appendix B: More details about the analytical model
The goal of this section is to present more details of the analytical model we introduced to capture the main behaviorof the metapopulation Schelling model presented in the main text.Let us denote by n Ai ( t ), respectively n Bi ( t ), the number of agents of type A , respectively B , at node i –th attime t . Each node is also characterized by a number of vacancies n Ei ( t ) at time t , however for all t and all i onehas n Ai ( t ) + n Bi ( t ) + n Ei ( t ) = L and moreover the number of agents of each kind is a preserved quantity, i.e. thesystem is closed. The state of the system is thus completely characterized by knowledge of ( (cid:126)n A ( t ) , (cid:126)n B ( t )), being (cid:126)n A ( t ) = ( n A ( t ) , . . . , n AN ( t )) and (cid:126)n B ( t ) = ( n B ( t ) , . . . , n BN ( t )).An agent A in node i –th is unhappy - or her fitness is low - if the fraction of B in the same node, that is n Bi / ( n Ai + n Bi ),is larger than a given tolerance threshold (cid:15) ∈ (0 , A is unhappy at i if n Bi n Ai + n Bi ≥ (cid:15) , (B1)and similarly for B .We assume that once an agent decides to move because unhappy, she will move to an uniformly randomly chosennew node, provide there is enough space there (that we named weak liquid Schelling model ). The model is intrinsicallystochastic and hence it can be described by the probability to be at time t in state ( (cid:126)n A , (cid:126)n B ), that is P ( (cid:126)n A , (cid:126)n B , t ).The evolution of such probability can be obtained using the master equation: P ( (cid:126)n A , (cid:126)n B , t + 1) = P ( (cid:126)n A , (cid:126)n B , t ) + (cid:88) ( (cid:126)n A (cid:48) ,(cid:126)n B (cid:48) ) (cid:104) T ( (cid:126)n A , (cid:126)n B | (cid:126)n A (cid:48) , (cid:126)n B (cid:48) ) P ( (cid:126)n A (cid:48) , (cid:126)n B (cid:48) , t ) − T ( (cid:126)n A (cid:48) , (cid:126)n B (cid:48) | (cid:126)n A , (cid:126)n B ) P ( (cid:126)n A , (cid:126)n B , t ) (cid:105) , (B2)1being T ( (cid:126)n A (cid:48) , (cid:126)n B (cid:48) | (cid:126)n A , (cid:126)n B ) the Transition probability to pass from state ( (cid:126)n A , (cid:126)n B ) to the new compatible one, i.e. thesystem can pass from the former to the latter, ( (cid:126)n A (cid:48) , (cid:126)n B (cid:48) ). Non compatible states cannot be linked and thus we willset T ( (cid:126)n A (cid:48) , (cid:126)n B (cid:48) | (cid:126)n A , (cid:126)n B ) = 0.The non zero transition probabilities can be computed using the behavioral rules of the model: A moves from node i to node j : T ( n Ai − , n Aj + 1 | n Ai , n Aj ) = 1 N n Ai n Ai + n Bi Θ (cid:18) n Bi n Ai + n Bi − (cid:15) (cid:19) n Ej LA moves from node j to node i : T ( n Ai + 1 , n Aj − | n Ai , n Aj ) = 1 N n Aj n Aj + n Bj Θ (cid:32) n Bj n Aj + n Bj − (cid:15) (cid:33) n Ei LB moves from node i to node j : T ( n Bi − , n Bj + 1 | n Bi , n Bj ) = 1 N n Bi n Ai + n Bi Θ (cid:18) n Ai n Ai + n Bi − (cid:15) (cid:19) n Ej LB moves from node j to node i : T ( n Bi + 1 , n Bj − | n Bi , n Bj ) = 1 N n Bj n Aj + n Bj Θ (cid:32) n Aj n Aj + n Bj − (cid:15) (cid:33) n Ei L . (B3)Let us observe that to lighten the notation we wrote only the variables whose values change because of the transition.The function Θ( x ) is defined to be 1 if x > • N : probability to draw the i –th node with uniform random probability; • n Ai n Ai + n Bi : probability to draw one A agent among the n Ai present over the total node population n Ai + n Bi ; • Θ (cid:16) n Bi n Ai + n Bi − (cid:15) (cid:17) : probability the selected A agent is unhappy; • n Ej L : probability to select a vacancy in node j .The master equation is unmanageable but one can go one step further by computing the time evolution of relevantquantities, for instance the average number of agents A and B in every node i at time t , (cid:104) n Ai (cid:105) ( t ) = (cid:80) (cid:126)n A n Ai P ( (cid:126)n A , (cid:126)n B , t ): (cid:104) n Ai (cid:105) ( t + 1) − (cid:104) n Ai (cid:105) ( t ) = (cid:88) (cid:126)n A (cid:88) j (cid:54) = i (cid:104) n Ai T ( n Ai , n Aj | n Ai + 1 , n Aj − P ( n Ai + 1 , n Aj − , t )+ n Ai T ( n Ai , n Aj | n Ai − , n Aj + 1) P ( n Ai − , n Aj + 1 , t ) − n Ai T ( n Ai − , n Aj + 1 | n Ai , n Aj ) P ( n Ai , n Aj , t ) − n Ai T ( n Ai + 1 , n Aj − | n Ai , n Aj ) P ( n Ai , n Aj , t ) (cid:105) = − (cid:88) j (cid:54) = i (cid:104) T ( n Ai − , n Aj + 1 | n Ai , n Aj ) (cid:105) + (cid:88) j (cid:54) = i (cid:104) T ( n Ai + 1 , n Aj − | n Ai , n Aj ) (cid:105) . And similarly for B : (cid:104) n Bi (cid:105) ( t + 1) − (cid:104) n Bi (cid:105) ( t ) = − (cid:88) j (cid:54) = i (cid:104) T ( n Bi − , n Bj + 1 | n Bi , n Bj ) (cid:105) + (cid:88) j (cid:54) = i (cid:104) T ( n Bi + 1 , n Bj − | n Bi , n Bj ) (cid:105) . Using the expression for the transition probabilities (B3) and assuming correlations can be neglected, i.e. (cid:104) ( n Ai ) (cid:105) ∼(cid:104) n Ai (cid:105) , we get: (cid:104) n Ai (cid:105) ( t + 1) − (cid:104) n Ai (cid:105) ( t ) = − N (cid:88) j (cid:54) = i (cid:104) n Ai (cid:105)(cid:104) n Ai (cid:105) + (cid:104) n Bi (cid:105) Θ (cid:18) (cid:104) n Bi (cid:105)(cid:104) n Ai (cid:105) + (cid:104) n Bi (cid:105) − (cid:15) (cid:19) L − (cid:104) n Aj (cid:105) − (cid:104) n Bj (cid:105) L + 1 N (cid:88) j (cid:54) = i (cid:104) n Aj (cid:105)(cid:104) n Aj (cid:105) + (cid:104) n Bj (cid:105) Θ (cid:32) (cid:104) n Bj (cid:105)(cid:104) n Aj (cid:105) + (cid:104) n Bj (cid:105) − (cid:15) (cid:33) L − (cid:104) n Ai (cid:105) − (cid:104) n Bi (cid:105) L . (B4)and (cid:104) n Bi (cid:105) ( t + 1) − (cid:104) n Bi (cid:105) ( t ) = − N (cid:88) j (cid:54) = i (cid:104) n Bi (cid:105)(cid:104) n Ai (cid:105) + (cid:104) n Bi (cid:105) Θ (cid:18) (cid:104) n Ai (cid:105)(cid:104) n Ai (cid:105) + (cid:104) n Bi (cid:105) − (cid:15) (cid:19) L − (cid:104) n Aj (cid:105) − (cid:104) n Bj (cid:105) L + 1 N (cid:88) j (cid:54) = i (cid:104) n Bj (cid:105)(cid:104) n Aj (cid:105) + (cid:104) n Bj (cid:105) Θ (cid:32) (cid:104) n Aj (cid:105)(cid:104) n Aj (cid:105) + (cid:104) n Bj (cid:105) − (cid:15) (cid:33) L − (cid:104) n Ai (cid:105) − (cid:104) n Bi (cid:105) L . (B5)2Observe that the right hand sides of (B4) and (B5) remain unchanged if we allow both sums to run over all nodeindexes, namely include also j = i .To go one step further, let us define the fraction of A and B in each node, that is α i = (cid:104) n Ai (cid:105) L and β i = (cid:104) n Bi (cid:105) L , (B6)then we can rewrite the previous equations (B4) and (B5) in terms of α i and β i . Finally rescaling time by s = t/L ,dividing the equations for α i and β i by 1 /L and passing to the limit L → + ∞ we get: dα i ( s ) ds = lim L → + ∞ α i ( s + 1 /L ) − α i ( s )1 /L (B7)= − N α i α i + β i Θ (cid:18) β i α i + β i − (cid:15) (cid:19) (cid:88) j (1 − α j − β j ) + 1 − α i − β i N (cid:88) j α j α j + β j Θ (cid:18) β j α j + β j − (cid:15) (cid:19) , and dβ i ( s ) ds = lim L → + ∞ β i ( s + 1 /L ) − β i ( s )1 /L (B8)= − N β i α i + β i Θ (cid:18) α i α i + β i − (cid:15) (cid:19) (cid:88) j (1 − α j − β j ) + 1 − α i − β i N (cid:88) j β j α j + β j Θ (cid:18) α j α j + β j − (cid:15) (cid:19) . Remark B.1 (Some preserved quantities) . Let us observe that the model correctly preserves the total fractions of A and B agents in time and thus the total vacancies (cid:80) i (1 − α i ( s ) − β i ( s )) = (cid:80) i (1 − α i (0) − β i (0)) = ρN , where ρ isthe emptiness defined previously, i.e. the total fraction of vacancies.To prove this statement is enough to take the time derivative of (cid:80) i (1 − α i ( s ) − β i ( s )) and observe that usingEqs. (B7) and (B8) one gets: dds (cid:88) i (1 − α i ( s ) − β i ( s )) = 0 . One can similarly prove the statement about the total fraction of A and B agents So in conclusion the system is ruled by the following system of differential equations: dα i ( s ) ds = − ρ α i α i + β i Θ (cid:16) β i α i + β i − (cid:15) (cid:17) + (1 − α i − β i ) N (cid:80) j α j α j + β j Θ (cid:16) β j α j + β j − (cid:15) (cid:17) dβ i ( s ) ds = − ρ β i α i + β i Θ (cid:16) α i α i + β i − (cid:15) (cid:17) + (1 − α i − β i ) N (cid:80) j β j α j + β j Θ (cid:16) α j α j + β j − (cid:15) (cid:17) . (B9) Appendix C: The average magnetization
To better analyze the system and in particular be able to describe the dependence of the average magnetization on (cid:15) , we introduce a new set of variables, the local emptiness , γ i = 1 − α i − β i , and the local difference of fractions of A and B , ζ i = α i − β i . The original coordinates can be obtained back using α i = (1 − γ i + ζ i ) / β i = (1 − γ i − ζ i ) / ζ i ∈ [ − ,
1] and γ i ∈ [0 ,
1] and the magnetization rewrites as (cid:104) µ (cid:105) = 1 N (cid:88) i | ζ i | − γ i . (C1)Let us introduce the functions F ( x ) := x Θ(1 − x − (cid:15) ) + (1 − x )Θ( x − (cid:15) ) and G ( x ) := x Θ(1 − x − (cid:15) ) − (1 − x )Θ( x − (cid:15) ) , (C2) Let us observe that (cid:104) n Ai (cid:105) / ( (cid:104) n Ai (cid:105) + (cid:104) n Bi (cid:105) ) − (cid:15) is positive if and only if α i / ( α i + β i ) − (cid:15) >
0, and similarly for the other term. α i α i + β i = 12 + ζ i − γ i ) and β i α i + β i = 12 − ζ i − γ i ) , we can rewrite Eq. (B9) as follows dγ i ( s ) ds = ρF (cid:16) + ζ i − γ i ) (cid:17) − γ i N (cid:80) j F (cid:16) + ζ j − γ j ) (cid:17) dζ i ( s ) ds = − ρG (cid:16) + ζ i − γ i ) (cid:17) + γ i N (cid:80) j G (cid:16) + ζ j − γ j ) (cid:17) . (C3)To qualitatively study the solutions of the latter system we define four zones in the ( ζ, γ ) plane as follows (see Fig.3 main text): Z = { ( ζ, γ ) ∈ [ − , × [0 ,
1] : 12 − ζ i − γ i ) − (cid:15) ≥ ζ i − γ i ) − (cid:15) < } Z = { ( ζ, γ ) ∈ [ − , × [0 ,
1] : 12 − ζ i − γ i ) − (cid:15) ≥ ζ i − γ i ) − (cid:15) ≥ } Z = { ( ζ, γ ) ∈ [ − , × [0 ,
1] : 12 − ζ i − γ i ) − (cid:15) < ζ i − γ i ) − (cid:15) ≥ } Z = { ( ζ, γ ) ∈ [ − , × [0 ,
1] : 12 − ζ i − γ i ) − (cid:15) < ζ i − γ i ) − (cid:15) < } . Let us observe that if (cid:15) < / Z is empty and if (cid:15) > / Z is empty.The initial condition of the Shelling model presented in the main text translate into ( ζ i , γ i ) distributed close to(0 , ρ ). We will thus be interested in the behavior of the solutions of Eq. (C3) in the Z zone if (cid:15) < / Z zone if (cid:15) > / (cid:15) < /
2, then one can easily prove that for all ( ζ, γ ) ∈ Z one has F (cid:16) + ζ − γ ) (cid:17) = 1 and G (cid:16) + ζ − γ ) (cid:17) = ζ − γ , so assuming that for all i one has ( ζ i (0) , γ i (0)) ∈ Z , then: (cid:40) dγ i ( s ) ds = ρ − γ idζ i ( s ) ds = − ρ ζ i − γ i + γ i N (cid:80) j ζ j − γ j , (C4)as long as ( ζ i ( t ) , γ i ( t )) will not leave Z . Assume for a while this statement, then the first equation can be straight-forwardly solved to give γ i ( t ) = ρ + e − t ( γ i (0) − ρ ), that is for all i , γ i ( t ) → ρ when t → ∞ .The solution for the second is more cumbersome, but one prove the existence of a lower (upper) solution ζ − i ( t )( ζ + i ( t )) such that ζ − i ( t ) < ζ i ( t ) < ζ + i ( t ) where: ζ ± i ( t ) = (cid:18) e − t (1 − γ i (0))1 − ρ − e − t ( γ i (0) − ρ ) (cid:19) ρ − ρ (cid:34) ζ i (0) ∓ (1 − (cid:15) ) (cid:18) γ i (0) − − γ i (0) (cid:19) ρ − ρ (cid:90) ρ + e − t ( γ i (0) − ρ ) γ i (0) x (1 − x ) ρ − ρ ( x − ρ ) ρ − ρ dx (cid:35) , (C5)from which one can get ζ i ( t ) → i when t → , ρ ) belongs to Z we have a posteriori proved that ( ζ i ( t ) , γ i ( t )) will never leave Z . A simplerbut also weaker statement can be obtained by observing that the equilibrium ( ζ i , γ i ) = (0 , ρ ) for all i is stable beingits eigenvalues − ρ/ (1 − ρ ) and −
1, each one with multiplicity N , thus there exists a neighborhood of (0 , ρ ) such thatall orbits whose initial conditions are inside never leave it. Using the definition of (cid:104) µ (cid:105) given by Eq. (C1), we have thusproven that for (cid:15) < . (cid:15) > /
2, then one can easily prove that for all ( ζ, γ ) ∈ Z one has F (cid:16) + ζ − γ ) (cid:17) = G (cid:16) + ζ − γ ) (cid:17) =0, so assuming that for all i one has ( ζ i (0) , γ i (0)) ∈ Z , then: (cid:40) dγ i ( s ) ds = 0 dζ i ( s ) ds = 0 , (C6)and thus ( ζ i ( t ) , γ i ( t )) will not evolve and thus remain in Z . However fixed the size of the domain centred at( ζ, γ ) = (0 , ρ ), where initial conditions are taken, the above assumption cannot be satisfied if (cid:15) > . Z shrinks to zero as (cid:15) → . (cid:104) µ (cid:105) also for (cid:15) > . F and G vanish for ( ζ, γ ) ∈ Z , thus Eq. (C3) can be rewritten as: dγ i ( s ) ds = ρF (cid:16) + ζ i − γ i ) (cid:17) − γ i N (cid:104)(cid:80) j ∈ Z F (cid:16) + ζ j − γ j ) (cid:17) + (cid:80) j ∈ Z F (cid:16) + ζ j − γ j ) (cid:17)(cid:105) dζ i ( s ) ds = − ρG (cid:16) + ζ i − γ i ) (cid:17) + γ i N (cid:104)(cid:80) j ∈ Z G (cid:16) + ζ j − γ j ) (cid:17) + (cid:80) j ∈ Z G (cid:16) + ζ j − γ j ) (cid:17)(cid:105) , (C7)where we used the shortened notation j ∈ Z k , to mean ( ζ j , γ j ) ∈ Z k and where the zones Z and Z have been definedin the Fig. 3 of the main text.The initialization of the metapopulation model, translates into original variables ( α i , β i ) initially uniformly randomlydistributed in (1 − ρ ) / δU [ − , δ . Such domain is distorted passing to the new variables( zeta i , γ i ), in particular it is translated into (0 , ρ ), rotated by 45 ◦ and expanded by a factor √ F (0) = G (0) = F (1) = G (1) = 0 and thus Eqs. (C7) admit as equilibrium points ζ i = ± (1 − γ i ).In conclusion, initial conditions inside Z will remain there while orbits originated from initial conditions in Z and Z will converge somewhere onto the straight lines ζ = ± (1 − γ ). In the left panel of Fig. 8 we report 5000 initialconditions built using the above described procedure for δ = 0 .
02, points in Z are marked in green, points in Z inblue and points in Z in red. The Eqs. (C7) are then numerically solved and the asymptotic positions are drawn inthe right panel of Fig. 8, giving to any points the color corresponding to the zone from where is started. One canclearly observe that all blue points remain in the Z zone, while green (respectively red) points converge to ζ = γ − ζ = 1 − γ ). ζ -0.15 -0.1 -0.05 0 0.05 0.1 0.15 γ ζ -0.15 -0.1 -0.05 0 0.05 0.1 0.15 γ FIG. 8: Dynamics for (cid:15) > .
5. Left panel: 5000 initial conditions are uniformly drawn ( α i , β i ) ∈ (1 − ρ ) / δU [ − ,
1] andthen transformed into the new variables ( ζ i , γ i ); points in Z are colored in green, points in Z in blue and points in Z in red.Right panel: asymptotic configuration of the orbits corresponding to the initial conditions given in the left panel, each point isdrawn using its initial color. The black solid line is the straight line ζ = (2 (cid:15) − − γ ), the dot-dashed line is the straight line ζ = (1 − (cid:15) )(1 − γ ) and the dashed lines (right panel) are the straight lines ζ = ± (1 − γ ). We are now able to provide a good approximation of the average magnetization. First of all let us rewrite Eq. (C1)as: (cid:104) µ (cid:105) = 1 N (cid:32) (cid:88) i ∈ Z | ζ i | − γ i + (cid:88) i ∈ Z ∪ Z | ζ i | − γ i (cid:33) = N N N (cid:88) i ∈ Z | ζ i | − γ i + N N N (cid:88) i ∈ Z ∪ Z | ζ i | − γ i , where N is the number of points in the Z zone and N the total number of points in the Z and Z zones. Becausepoints in Z and Z converge to ζ = ± (1 − γ ) the right most term in the previous equation is trivially equal to N /N .To compute the contribution arising from points in Z is more cumbersome, let us hereby stress the main ideas. Onecan easily show that if (cid:15) > . δ/ (1 − ρ ) then the lines ζ = (2 (cid:15) − − γ ) and ζ = (1 − (cid:15) )(1 − γ ) do not intersect thediamond like domain (see left panel Fig. 9) and thus all points are initially in the Z zones, this means that N = N and N = 0 so the magnetization has only a contribution from the Z zone. On the contrary if (cid:15) < . δ/ (1 − ρ )5 ζγ ε -1+2 ε ζγ ε -1+2 ε QP1 ε>1/2 + δ/(1−ρ) ε<1/2 + δ/(1−ρ)
P’Q’ T= (0,ρ+2δ) T= (0,ρ+2δ) R= (0,ρ−2δ) R= (0,ρ−2δ) FIG. 9: The geometries used to compute an approximation to the average magnetization in the case (cid:15) > . some initial conditions are in the Z and Z zones (see right panel Fig. 9) and thus the magnetization has bothcontributions.In the case (cid:15) > . δ/ (1 − ρ ) one can estimate the contribution to the magnetization as: µ int, = 1 N (cid:88) i ∈ Z | ζ i | − γ i ∼ δ (cid:90) δ dζ ζ (cid:90) − ζ + ρ +2 δζ + ρ − δ dγ − γ = 14 δ (cid:90) δ dζ ζ log 1 − ζ − ρ + 2 δ ζ − ρ − δ , (C8)being 8 δ the measure of the blue diamond domain in the left panel of Fig. 9. While for (cid:15) < . δ/ (1 − ρ ) one canfind the following estimation µ int, = 1 N (cid:88) i ∈ Z | ζ i | − γ i ∼ P (cid:90) ( ζ P + ζ Q ) / dζ ζ (cid:90) − ζ + ρ +2 δζ + ρ − δ dγ − γ = 2 P (cid:90) ( ζ P + ζ Q ) / dζ ζ log 1 − ζ − ρ + 2 δ ζ − ρ − δ , (C9)where P is the measure of the polygon T P QRQ (cid:48) P (cid:48) (in blue in the right hand side of Fig. 9). Let us observe thatin the last integral we make the approximation that points P and Q can be replaced by an average point whose ζ coordinate is the average of the ones for P and Q ; this is a minor assumption that helps to compute the integral andwill not influence the final result as show below. Now both integrals can be exactly computed.To get the final estimate for (cid:104) µ (cid:105) we need to compute N /N and N /N in the case (cid:15) < . δ/ (1 − ρ ), the othercase being trivial and already considered. Assuming the number of points sufficiently large we can affirm that suchfractions are well approximated by the ratio of the corresponding polygons, that is N N ∼ P δ and N N ∼ − P δ . In conclusion we can approximate the average magnetization for all (cid:15) > . (cid:104) µ (cid:105) approx = P δ µ int,k + (1 − P δ ) , (C10)where µ int, is given by Eq. (C8) valid for (cid:15) > . δ/ (1 − ρ ), and µ int, by Eq. (C9) for the case (cid:15) < . δ/ (1 − ρ ).In Fig. 10 we compare the average magnetization obtained using its very first definition Eq. (C1) and the numericalintegration of the system (C3), with the approximation given by Eq. (C10) for several values of δ and ρ = 0 . ǫ < µ > FIG. 10: Averaged magnetization (cid:104) µ (cid:105) and the approximation (cid:104) µ (cid:105) approx as a function of (cid:15) > .
5. Symbols correspond tonumerical integrations of the system (C3) and the use of the definition for the magnetization, lines to the calculation of theapproximation Eq. (C10). Blue circles correspond to δ = 0 .
02, black squares to δ = 0 .
01 and red diamonds to δ = 0 . ρ = 0 ..