Self-Organized Societies: On the Sakoda Model of Social Interactions

We characterize the behavior and the social structures appearing from a model of general social interaction proposed by Sakoda. The model consists of two interacting populations in a two-dimensional periodic lattice with empty sites. It contemplates a set of simple rules that combine attitudes, ranges of interactions, and movement decisions. We analyze the evolution of the 45 different interaction rules via a Potts-like energy function which drives the system irreversibly to an equilibrium or a steady state.We discuss the robustness of the social structures, dynamical behaviors, and the existence of spatial long range order in terms of the social interactions and the equilibrium energy.


Introduction
Throughout the history, the society has self-organized in a myriad of social structures and behaviors which appears as a response to the attitudes, decisions, and expectations among individuals, that is to say, from local (usually simple) rules to global behavior [1].This capacity appears in animals [2], insects [3,4], among others, and it flourishes in human beings since the early life [5].The basic underlying mechanism of this self-organization lies on the evaluation of the expectations resulting from the attitudes among individuals inducing, or not, a mobility decision.
Interestingly, despite the complex nature of social interactions, the social behavior shares many common features with a variety of physical systems.Segregation, inclusion, and aggregation are examples of a collective order arising from simple and local individual rules based on attitudes, decisions, and expectations among individuals [6][7][8][9].
In the late forties, Sakoda introduced in his Ph.D. thesis [10] a first general discrete dynamical model for social interactions and published later [11], at the same time of the well-known Schelling's social segregation model [12,13], which is, as we will show, only a special class already included in the original Sakoda model.The particular case of the Schelling segregation model and its variations have been studied widely through agent based models [9,[14][15][16][17][18][19][20][21][22][23].However, the Sakoda dynamics and its underlying richness, which explain other different social phenomena far from segregation, are mostly unknown; thus they have not been, to the best of our knowledge, studied previously in a deeper manner.
The original Sakoda model roughly consists of social interactions among two groups of individuals evolving in a network according to specific attitudes of attraction, repulsion, and neutrality.An individual evaluates its social expectative at all possible available locations, preferring those near individuals associated with attractive (positive) attitudes and avoiding locations near individuals associated with repulsive (negative) ones.This procedure is repeated randomly among all possible individuals; henceforth Sakoda's algorithm is iterated recursively driving the system, under some conditions, to a well-organized spatial pattern.
The dynamics is quite rich because of three aspects: (1) the large number of combinations of the possible attitudes, (2) the effect of the separation distance among individuals (i.e., the interaction could be of short or long range), and (3) the 2 Complexity mobility of the individuals (the individuals may move into its own neighborhood or they may migrate far away).This produces a myriad of possible patterns of the individuals' arrangements in the space (social structures) to be analyzed.Some of them were already recognized by Sakoda in his early work [11].
In the present paper we discuss and complement the Sakoda original classification of social structures.More precisely, we characterize all independent individual interactions (45 different interactions) in the case of short range interaction and long range movement.We introduce a Potts energy [24,25], which has been the natural -state generalization of the Ising model.As the former, the Potts model also displays a phase transition in statistical physics.In the context of the Sakoda model, the Potts energy happens to be a Lyapunov function in the case of symmetric interactions that governs the future dynamics and quantifies different attractors in terms of the parameters of the problem.This Potts-like energy turns out to be a common principle for all possible interactions and describes the evolution of different social interacting systems.Finally, we discuss the role of empty sites and the existence of long range order in the social structures.
The present paper is organized as follows.Section 2 presents the Sakoda model and the numerical scheme.Section 3.1 shows the energy principle that mandates the dynamical evolution of the system.Next, Section 4 characterizes the various social structures arising from different interactions possibilities.A quantitative energy-based study is done in Section 5.1.The role of vacancies is studied in Section 5.2 and the existence of long range order is discussed in Section 5. 3. Finally, we conclude in Section 6 and complement the analytics with a series of appendices.

The Basic Model
The Sakoda model consists in a regular and periodic lattice with  nodes in a two-dimensional space.Each node  is associated with a discrete variable   and may take the values −1, 0, +1.States ±1 denote an occupancy by an individual of one of the two types (+ or −); the state 0 denotes an empty place or a vacancy.Only when a node  is empty (  = 0) may an individual (+ or −) occupy the free spot and reassign its value with its current one.
The social interaction is mainly characterized by the possible attitudes among individuals which are summarized in a 2 × 2 matrix.We call hereafter the "-matrix," which for the sake of space we will symbolize in a vectorized form:  = ( 11 ,  12 ,  21 ,  22 ).The entries   take three possible values: {+1, 0, −1}.These indicate, respectively, an attractive (+1), a neutral (0), or a repulsive (−1) attitude from members of the type  to the individuals of the type .There exist 3 4 = 81 different possible cases which are reduced, because of symmetry, to 45 independent social structures.In Appendix A, we discuss the symmetries of the "-matrices" and its properties in a deeper manner.
Next, we model the preferences of one individual over a particular place in the lattice based on the spatial locations and the attitudes of other individuals.Following the original ideas of Sakoda, we propose a function that quantifies the social expectations of the individual : where   ≥ 0 denotes a symmetric (  =   ) interaction strength.We avoid self-interactions taking   = 0. On the other hand, the sign of the interaction is given by (do not confuse the variable   that takes the values   = 0, ±1 with the entries of the matrix   that take the same values   = 0, ±1): In general, this expression is not symmetric; that is,   (  ,   ) ̸ =   (  ,   ).We underline that though   is symmetric, the full interaction     (  ,   ) is not necessarily symmetric.The coefficients   may include short and long range interactions.The mobility of an individual could be also of long or short movement.An individual located at the node  would move towards an empty node , if   (  ) <   (  ).
Here   (  ) represents the potential value (1) evaluated at the empty node, , occupied by the individual   .As a general rule, the individual chooses the place  everywhere in the lattice that makes the social expectation function   (  ) the highest.In case of degeneracy, that is, if there are  nodes,  1 ,  2 , . . .  , such that   1 (  ) =   2 (  ) = ⋅ ⋅ ⋅ =    (  ), one of these is selected randomly with the same probability.In a short movement case, the mobility of the individual is restricted to its Moore neighborhood (the eight closest neighbors).
To summarize, the movement of individuals occurs in the following manner: during an iteration step, an individual is selected randomly.Then, it evaluates the potential function at every available empty site in its range of movement for a potential new location.Next, the individual selects the site with the highest value of the potential expectation (1) and then it moves to this site.If there is not a site that improves the potential function, then the individual stays in place.The algorithm is iterated until the system reaches a fixed point; otherwise, it runs forever.
We point out that the rule preserves the initial number of individuals of all types.That is, if  + ,  − , and   are the number of individuals of the types +, −, and the number of vacancies, respectively, they keep their initial value during the evolution.We define the fractions by  ± =  ± / and   =   / (naturally  + +  − +   = 1).Hereafter, as in Sakoda original paper [11], we limit our study to the case of  + =  − = (1 −   )/2.

Numerics.
For the following, we restrict our central results only to the case of short range interactions and long range movements: an individual evaluates its preferences considering the Moore neighborhood; that is,   = 1, if  is in the 8 closest neighbors of the site , and   = 0 elsewhere.Then, the individual may look forward the highest expectation at any empty node in the lattice.We discuss the cases of short and long range interactions/mobility at the end of the paper and in Appendix F.
The numerical simulations for all 45 possible -matrices were run in a two-dimensional periodic lattice of  = 128 × 128.The evolution usually runs for 2 × 10 5 time steps, which happens to be sufficient to reach equilibrium or a steady state.We have also simulated the case of  = 64 × 64 without significant dependence on the system size.
The initial condition consists of a random distribution of states characterized by a fraction   of vacancy states.For all of our simulations, we have checked the robustness of the social structures considering different random initial configurations.We focus our attention on five levels (  = 0.01, 0.25, 0.5, 0.75 & 0.98) of different initial concentration of vacancies as a good representative sample (see Section 5.2).Other initial distributions may be studied using the same frame and tools developed in the present paper.From the simulations, we notice the similarities of our results with the physics of surface tension [26].Indeed all possible situations among the three phases follow after an interface equilibria argument (something that will be clear in Section 3).For instance, Figure 1(a) shows a situation of a miscible +1 and −1 phases, but both are nonmiscible with respect to white.Similarly, Figure 1(b) presents a miscible phase between +1 and 0 (white), but both are nonmiscible with −1 phase.Figures 1(c) and 1(d) present cases with the three nonmiscible phases (+1, −1, 0), but the interface equilibria are different.While in Figure 1(c) there are interfaces between +1/−1, +1/0, and −1/0; in the case of Figure 1(d), +1/−1 Complexity interface does not exist and it is replaced by a monolayer of white vacancies (excepting the case of a very rare number of vacancies, see Figure 1(d)(I) for   = 0.01).For a given interaction, a robust behavior is clearly noticed for different vacancy fractions, but the patterns may differ substantially from one case to another because of the role of the empty sites.We shall discuss more in detail this point later on in Section 5.2.
Finally, we have checked the robustness of the phenomena for the case of a nondeterministic Sakoda instance, through considering the rule: An individual located at node  would move towards an empty node , with a probability   (  ) = (1/) (  (  )−  (  )) .Here  = 1 + ∑   (  (  )−  (  )) .As expected, in the large  limit one recovers the deterministic rule; on the contrary, in the limit  → 0 the individuals move randomly to any available empty site.In this case the system behavior is highly fluctuating and it does not selforganize into any spatial structure.Therefore, it is expected that, for a given -matrix, a critical parameter   exists, which controls the existence or not of a social structure.

An Energy Theorem for the Dynamics of the Sakoda Model.
The evolution of the Sakoda model is ruled by an energy principle that we formulate in this section.After [27], the following theorem is shown.Theorem 1.If the -matrix is symmetric, then the Potts-like energy function [24,25], of a given configuration of the network  = { Rewriting the r.h.s. of (4) in terms of the expressions   and   of ( 1) and ( 5), respectively, we obtain Δ =   (  ) −   (  ) ≤ 0. We emphasize that the equality holds in the case of   (  ) =   (  ); in this case no exchange is carried out.In conclusion, if  is a symmetric matrix, then Δ ≤ 0; that is, energy (3) is a nonincreasing quantity after a movement, concluding the proof.We emphasize that the proof is completely general; hence, it does not depend neither on the range of the interactions and the range of motion.
If the -matrix is not symmetric, energy (3) is not necessarily a Lyapunov function; thus the energy may increase or decrease indistinctly.However, as we shall see next, there is an energy decreasing principle for both cases, symmetric and nonsymmetric interactions.This principle is understood in a statistical sense: for a random configuration, in average, it is more probable to decrease than to increase the energy after a movement [23].
This principle is understood in a statistical sense: for a random configuration, in average, it is more probable to decrease than to increase the energy after a movement [23].The observation of an energy decreasing principle for nonsymmetric interactions may indicate behind this principle a "free energy."In [21] we concluded that this free energy is  =  −  +  + −  −  − −     , where  ± and   are the Lagrange multipliers associated with the constant number of individuals of each type and the constant number of vacancies.Hence, the macrobehavior minimizes energy with the corresponding restrictions (see next section).
Finally, notice that the energy (3) may be written in terms of geometrical parameters [16,20,21].For instance, in the case considered here of short range interaction, the energy reads where ℓ  =   − (1/3)  is an effective interface length which is directly related to the number of edges,   , and corners,   , among the interface created by the states  and  (see Appendix D for a derivation of ( 6)).Finally, the term −4( 11  + +  22  − ) represents the bulk energy contribution.However, because  ± =  ± are constant during the evolution, the bulk energy term does not change; therefore variations of the energy are produced only by interface length variations, if these exist.As we will see later, the parameter  plays an important role in the social structure characterization.
(Color online) temporal evolution of the energy (i) and the effective perimeters ℓ  (ii) for the same cases as in Figure 1 row (II), that is, for   = 1/2.The curves are labeled with an (a) for the snapshot (a II) of Figure 1, with a (b) for (b II) and so forth.(i) In this plot one notices clearly energy decreasing tendency in time.Here the curves (b) and (c) correspond to a nonsymmetric -matrix; hence its energy may decrease or increase indistinctly.Some increments of energy are labeled with * .(ii) The perimeters ℓ  are explicitly labeled in the plot with +−, +0, or −0.It is observed that the qualitative evolution agrees with the energy behavior.These special cases ran for 1.8 × 10 6 time steps.

Energy Behavior in the Sakoda Model.
To illustrate the energy principle, Figure 2(i) plots the temporal evolution of the energy (3) for the same cases of Figure 1.The selection includes both symmetric and nonsymmetric interactions.It is easy to check that this energy function is a Lyapunov function for the former cases, while for the later it is not (see Figure 2(i)).Nevertheless, for all nonsymmetric interactions, the energy tends to a certain value indicating the existence of an equilibrium state (see curves (a) and (d) in Figure 2).This case suggests the existence of some kind of "statistical attractor," something already noticed in Schelling's model [23].Despite the complexity involved in the Sakoda model, the long-time dynamics reaches a statistical attractor with a well-defined average energy and social structure (see Section 5.1).
Quantitatively, the energy expansion (6) allows us to distinguish between the segregation and aggregation patterns.Indeed, the energy contributions , and E −0 = (3/2) 22 ℓ −0 correspond, respectively, to the interface energies among the two populations + and −, +, and 0 state, and − and the vacancies.In complete analogy with the physics of surface tension [26], the coefficients ,  11 , and  22 would be a kind of "social tension" among the interfaces.From a physical point of view, an interface is nonmiscible if the respective social tension, namely, ,  11 , or  22 , is positive.Hence, because the system minimizes the energy, it minimizes the corresponding perimeter; otherwise, the interface would be miscible.Thereby the interface increases its length mixing the respective states.
The previous arguments allow us to give an explanation about why segregation is more efficient in the case of nonsymmetric interactions; for example, the segregation is stronger in Figure 1(c)(II), than in the symmetric case Figure 1

(d)(II).
In the nonsymmetric interaction case, the energy is not mandated strictly to decrease (see Figure 2(a)).Therefore, the system may evolve forever, searching in average the lowest energy as possible [23].On the other hand, in the symmetric interaction case, the system reaches a frustrated state with a larger energy.
As a consequence of this, the whole patterns of Figure 1 may be understood according to (6).In Figure 1(a), that is,  = (0, 1, 1, 0) and  = −2, the perimeter energy reads −3ℓ +− .In this case, the energy decreases if the interface perimeter between −1 and +1 populations increases (see ℓ +− blue curve in Figure 2(b)), thus creating a labyrinth pattern in analogy with patterns of antiferromagnetic materials.
In the case of Figure 1(b), that is,  = (1, 1, −1, −1) and  = 0, the energy term is simplified to (3/2)(ℓ +0 − ℓ −0 ), being this a minima if, simultaneously, ℓ +0 is maximized (see the upper red curve in Figure 2(b)) and ℓ −0 is minimized (see the ℓ +0 red curve in Figure 2(b)).In the former case, +1 individuals arrange in stripe patterns maximizing the perimeter with the white vacancies, while in the latter case, −1 individuals make clusters surrounded by a white layer of 0. Cases of Figure 1(c),  = (1, 1, −1, 1), that is,  = −2, and Figure 1(d),  = (1, −1, −1, 1) ( = 4), produce segregation patterns.For these cases, the perimeter energy terms are similar: 1(c) and 1(d), respectively.Both cases tend to minimize all perimeters among interfaces (see the green and orange curves in Figure 2(b)), but a notorious difference appears: in Figure 1(d) there is no +1/−1 interface.As it can be seen in Figure 2 (curve ℓ +− (d)), the interface perimeter, ℓ +− , vanishes as the time goes by.This happens because forming this interface "costs" twice compared to forming the same interface of Figure 1(d); this produces that +1 and −1 regions are surrounded by a white layer which has no cost.

The Basic Social Structures
Under the same conditions of Figure 1, we have run numerical simulations for all 45 possible -matrices.The social Table 1: Summary information on the patterns of Figure 3.The table is indexed by  1 and  2 .Each box contains the precise S-matrix and the exponent  obtained from the corresponding Fourier spectrum (see Section 5.3).The cases labeled " = none" mean that the slope may not be defined.3 and in Table 1.
The phase diagram of Figure 3 allows us to classify the basic social behaviors characterized by well-defined social structures.In his original work, Sakoda [11] summarized this rich dynamics by eight relevant behaviors with a particular interpretation of the real life situation (we keep original Sakoda's terminology): crossroads, mutual suspicion, segregation, social workers, boys and girls, couples, and husband and wife.Based on the symmetries of the diagram, we suggest in the present work three new social behaviors: Inclusion, Repulsion, and Exclusion, which complete the whole picture.Naturally, the Sakoda model behavior is only simplified view of real social behaviors that may be absent in the current model.
We underline that all patterns as well as its evolution may be understood quantitatively in terms of the energy minimization principle introduced in the present paper.In particular, as expected, the parameter  =  11 + 22 −( 12 + 21 ) splits roughly the segregation-like and the aggregation-like pattern.
The description of the social behaviors presented in Figure 3 and its respective social structures are as follows: (S1) Segregation typifies attraction among individuals of the same group and repulsion among members of the other group ( = 4).In this case, the system minimizes the interface perimeter inducing segregation (see Figure 1(d)).
(S6) Husband and wife is the case where the members of one group feel attraction by everybody, while individuals of the opposite group feel attraction for the former group and repulsion for individuals of the same group.In this case  = (1, 1, 1, −1).
The two previous cases correspond to  = 0, which characterizes a "mixed phase"; indeed, though the former, Inclusion, maximizes the interface perimeter ℓ +− , the later, Individualism, minimizes it.
(S10) Exclusion is a case where the individuals of one group want to be accepted by the opposite group, but the later do not accept the former group.Its -matrix is  = (−1, 1, −1, −1).

Quantitative Description of the Social Structures
5.1.Phase Diagram in the - Plane.All patterns may be understood quantitatively in terms of the energy minimization principle presented throughout this paper.Therefore, the energy seems to be a pertinent quantity for a characterization of the resulting social structures, at least concerning the social segregation [28].For all cases, the energy reaches a steady state allowing us to characterize the social structures as statistical attractors.
In the special case  ± = (1 −   )/2, the energy expression (6) reads ( ± =  ± ) From the numerical simulations, we conclude that the behavior of the corners and the interface length scale proportional also to (1 −   ), that is, the total energy scales, as / ∼ (1 −   )Φ 1 (,  1 ), where Φ 1 (,  1 ), is a function of  &  1 .The relationships among perimeters, corners, and the energy are more exhaustively developed in Appendix E. Figure 4 shows a set of 40 of the 45 attractors (the missing 5 social structures are overlapped) classified in the - plane for   = 1/2.This figure shows that the energies are below an upper envelope, labeled by (i) in Figure 4.It may be noticed the existence is an absolute energy maximum for all configurations for  = (−1, −1, −1, −1) that correspond to  = 0 and  1 =  11 +  22 = −2 and an absolute minimum  = (1, 1, 1, 1) ( = 0,  1 = 2).The social structures are spanned all over the - plane below this envelope (i) and above the line (ii) for  < 0 and the horizontal line /4(1 −   ) = −1 for  > 0. The existence of two different behaviors for  < 0 (aggregation) and  > 0 (segregation) is a clear consequence of the existence of these two phases.
Figure 4 presents a more complete landscape of the social structures than Figure 3, showing that energy appears to be a good characterization of social structures.However, the parametrization includes the energy that itself is a function on the interaction parameters.The energy spreads successfully the large number of degeneracies of the restricted classification in planes  1 and  2 of Figure 3.However, the energy itself depends on the parameters of the problem ( 1 ,  2 ,  3 ) (see (A.3)); therefore, the (, ) characterization is visually attractive, but the ( 1 ,  2 ) characterization reveals out the real symmetries of the interactions and its consequences.
Lastly, we notice the existence of the 11th social structure, labeled by S11 and named as Neutrals.The social structures belong to Neutrals class where the individuals of different groups do not express neither rejection nor attraction among others.This behavior is represented by the interaction  = (0, 0, 0, 0).Similarly, we identify the 9th group of social behavior (G9) named by social neutrality.

The Role of the Vacancy Fraction on the Social Structures.
As it can be noticed in Figure 1, the final patterns depend on the empty spaces in the lattice,   .In particular, we observe that the number of representative social structures diminishes as both   → 0 and also in the limit   → 1.In the case of large number of vacancies, the individuals move almost freely through the space; as the vacancies diminish, finite space effects play a role.This is explained because of the smaller number of options that individuals have to move in the lattice.This effect modifies the basic interaction dictated, a priori, by (1).For instance, an individual that feel repulsion in free space would not necessarily follow a repulsive movement in a restricted environment (see Figures 1(c)(I) and 1(d)(I)).We investigated systematically the effect of   on the final patterns.
Figure 5 summarizes the phase diagram of some attractors in the  −   space.The reader may note that different vacancy concentrations,   , do affect dramatically the final attractors.Nevertheless, one notices that the lines  = −2 and  = 1 characterize two critical lines which separate clearly three domains.For −4 ≤  < −2, we observe a heterogeneous phase, without long range order.For −2 <  < 1, a neutral phase appears with coexistence of a disordered phase.Lastly, for 2 ≤  ≤ 4, a segregation phase with quasi-long range order appears.In particular, we notice that, in the special case of Schelling's model, that is,  = (1, −1, −1, 1), one observes that if we increase the number of total vacancies we recover the phenomenology already found in [20].
The quasi-long range order is readily quantified by the spatial correlation function which reads () = (1/) ∑     + .It may be computed directly in terms of the Fourier transform of the spatial pattern: Here   is the spatial pattern,  = √  is the linear size of a square lattice, and  is an integer ranging within −/2 + 1 ≤  ≤ /2.A straight-forward calculation relates the correlation function with the spectrum.In the limit ,  → ∞, we have thus yielding In practice, we compute directly the fast Fourier transform (FFT) of the spatial patterns and after performing an angular average, one gets the isotropic spectrum (in ): The spectrum allows us to distinguish quantitatively the existence of different social structures; moreover, a sudden change of these correlations may denote a phase transition from one regime to another.For instance, in the segregation regime that corresponds to a ferromagnetic case with welldefined domains, one expects that   ∼  − , where  is a critical exponent.This power law behavior is a characteristic of quasi-long range order spatial correlation.Figure 6 presents the spectra of the four patterns of Figure 1 for   = 1/2.One notices that in the cases of (b), (c), and (d) one gets quasi-long range order with the following critical exponents:  (b)  = 2.6; (c)  = 2.9; and (d)  = 3.6.Notice that (c) possesses an exponent close to  = 3 which corresponds to Porod's law, indicating that the observed phenomena is related to phase separation [29].Besides the slight differences among the spectrum exponents, cases (b), (c), and (d) are systems that possess quasi-long range order, implying that faraway individuals are correlated among them.On the other hand, in the case of Figure 1(a), the spectrum is almost flat in the long length scales (short wave number ) but it shows a peak at the scale of the discrete lattice, revealing a short wavelength modulation (the stripe pattern).Table 1 summarizes the exponents for all cases corresponding to Figure 3.

Discussion
We have reported the dynamics and the social structures of the, almost forgotten, social interaction model proposed by Sakoda in the 1940s.The model is based on two different groups of individuals moving through a background of empty sites.It includes different attitudes that individuals assume (attractive, neutral, or repulsive) towards their own group and towards the other group, resulting in 45 fundamental interactions.All of them were studied in the present paper.We observe that different attractive/repulsive/neutral social behavior was mainly characterized depending on the parameters  11 and  22 and  = ( 11 +  22 ) − ( 12 +  21 ).These parameters allow us to classify the pattern by its key features.The entire spectra of possible attractors may be studied in the frame of well-known physical analogies, like ferromagnetic (segregation) and antiferromagnetic (heterogeneous attraction) patterns, kinetics ordering, phase transitions, existence of quasi-long range order, and surface tension phenomena.In this sense, we propose a Potts-like energy that turns out to be a Lyapunov function in the case of symmetric interactions.This energy characterizes different attractors in terms of the parameters of the problem, namely, the initial fraction of vacancy places, the range of interaction, and the range of motion.
Surprisingly, the Sakoda model amalgamates its underlying behavior under an energy principle which, under the condition of symmetric interactions, mandates the system to equilibrium.This energy appears to be a common principle for all possible interactions because it decreases statistically in time, beyond the range of validity of symmetric interactions.
In conclusion, typical social structures and behaviors arise as a consequence of the interactions among individuals and the migration to a more attractive place.Depending on the intrinsic interactions we observe a variety of patterns (see Figures 3 and 4); some of them may be a manifestation of our societies.Despite the complexity of the interactions, expectative, and motions, an energy minimization principle surges as common law that rules the evolution of social structures.We suggest that the pertinence of the energy principle may deserve more attention from the social science community.We underline that the present model is only an approximative view of social behaviors.As an example, all individuals are considered of the same type, they are not distinguishable, and all of them are simply characterized by a single -matrix.Real societies are certainly more complex and harder to understand.The aim of this paper is to stimulate further quantitative description based on Sakoda's approach.0.1 2.9 We end this paper commenting on possible generalization of the Sakoda model studied here: (1) There are many variations of this model that may be considered.First, our results may be generalized to any complex network of interactions, more general graphs, and other space dimensions.A second possibility is to attribute a nonzero entry in the -matrix for the interactions among states and vacancies.In this case, the -matrix becomes a 3 × 3 matrix, and the total number of possible interactions would be quite large: 3 9 = 19683, but only 3903 of them are independent.Similarly, in the case of -states of the variable   , the -matrix becomes  ×  interaction matrix and the total possible number of different cases to be considered would become huge: 3  2 .In practice, in these cases, only a small sample of social structures could be studied.Therefore, an energy unified picture, as the one presented here, would be helpful.
(2) Though we focused our work on the local range of influence of individuals and the long range of movement, the dynamical behavior for different ranges of influence and different movement of the individuals should be considered.Both the interaction range and the range of movements influence macroscopically the formation of the social structure.As a general rule, the typical length scale of the pattern depends on the range of interaction and/or the scale of the movement.The largest pattern arises for the case of long range movement/long range interaction.In this sense, agglomeration sizes are produced by the available information that individuals have access: individuals may explore more empty sizes if they consider a long range movement; on the contrary, short movements allow individuals to explore just few options of empty sizes; hence small clusters are expected.All possible cases of long/short interaction and long/short movement will be discussed in a separate publication.
Figure 7 summarizes the attractors for the same matrices of Figure 1 but for the cases of long movement/long range interactions and short range movement/short range interactions and short range movement/long range interactions.As a general rule, the typical length scale of the

C. Remark on the Potts Energy (3)
The -function equation ( 2) may be written as (we acknowledge J. P. Nadal for this remark): Notice that in the second equality the first three terms are symmetric under the exchange  ↔ , but the last one is not.After (C.1), the energy equation (3) takes the form of a spin-1 Ising model [20]: Because this energy is a symmetric sum on indexes  and , the antisymmetric terms of (C.1) cancel out.

D. Derivation of the Energy-Perimeter Formula (6)
The  terms in expression (C.2) contribute if and only if   and   are both different from zero.Hence, the sums run over ∑ → ∑  ≡ ∑ ,/    ̸ =0 producing, approximately, an effective sum only over (1 −   ) sites.In the case of Moore neighborhood, the edges and corners formula stated in (6) may be derived after the following considerations: let us define regions D + , D − , and D 0 (which are not necessarily simply connected); the regions composed by the sates are +1, −1, and 0, respectively.These domains possess  + ,  − , and  0 individuals.Further, we define  +− ,  +0 , and  −0 as the interfaces among D + /D − , D + /D 0 , and D − / D 0 ; these interfaces consider the presence of different states in the same neighbors.The interfaces are geometrically expressed in terms of edges and corners of the sates +1 and −1.(See Figure 8 as a sketch of the interfaces.) The energy is nonzero only if   and   are in D + or D − .For the individuals in the bulk of D + (the individuals   = +1 and whose neighbors are all   = +1) the energy contribution becomes where   + differs from  + because of the boundary or interface terms.Rewrite this expression in terms of Similarly, for the individuals   = −1 and whose neighbors are all   = −1, the energy becomes At the interface, the contribution on a site, say +1, to the energy comes from the eventual ++, +−, and +0 interaction.Assume that, in its neighborhood,   has  + individuals +,  − individuals of sign −, and  0 individuals at 0 (naturally  + +  − +  0 = 8), Then, after (C.2), the energy at an interface +− becomes In the case of a linear interface +−, the number of neighbors +1 is  + = 5,  − = 3, and  0 = 0. To summarize, in the case of a linear interfaces one gets Finally, one notices that as an effect of the discrete nature of the lattice, an interface may be composed of vertical and horizontal edges (joined by a corner).Hence, the total perimeter ℓ  is formed by the contribution of edges   and corners   .A careful examination provides the relation for the perimeter ℓ  in terms of edges and corners.One concludes that ℓ  =   − (1/3)  ; hence, we readily get (6).

E. The Edge-Corners Relationships
The numerical simulations provide evidence of correlations among the number of edges and corners of an interface for all concentration vacancies   together.Figure 9(a) plots the edges ( +− ,  +0 ,  −0 ) versus the respective corners ( +− ,  +0 ,  −0 ) for the 45 Sakoda cases.This linear relation among the perimeter and the corners has only a statistical sense which follows from the obvious equivalence among the number of edges and the number of corners in a diagonal interface.Naturally there are patterns with a large/small perimeter but with a small/large number of corners which are out of this description; for example, the heterogeneous attraction,  = (−1, 1, 1, −1), and the repulsion  = (−1, −1, −1, −1) are far away from the trend line.

F. The Long Range Interaction Case
In this case, let us consider the interaction exchange coefficient   = V(  ) a function that depends explicitly on the spatial distance among sites  and , being   = |  −   |, the Euclidean distance.Therefore, (2) of the main text may be written as To avoid any problem with a divergent self-interaction, more precisely, case  =  in sum (F.1), assume that V(0) = 0 in our model.Though V(⋅) could be any kind of function, as a general rule, consider that V(⋅) is a decreasing function in its argument; that is, the shorter the distance is, the stronger the influence is.
In the present chapter, for the long range interaction, the potential function is proposed as follows: Originally, Sakoda reported his results considering V(  ) = 1/ 1/2  , because he argued that the weight contribution of faraway individuals should be quite strong; that is, the social interactions are of "extremely" long range.In general, because strong long range interactions do not evolve to an equilibrium (at least in short time) these cases are avoided and limited to "weak" long range interactions, say V() ∼ 1/ 2 in two-dimensional space.In this case it is expected that the robust behavior does not depend explicitly on the detail of the individual interaction, something already noticed by Schelling [1].
Figure 10 presents a comparison between the case on long range interaction/long range mobility in the particular case of V() = 1/ 2 and the Sakoda strongly long range interaction V() = 1/√.The simulations show that the V() = 1/√ case does not depend strongly on the size of the system that is confirmed comparing Figure 7 (long influence/short movement case) with Figure 10.

Figure 5 :
Figure 5: (Color online) examples of different attractors classified by  parameter for different values of   as indicated in the plot.Numerical simulations run in a two-dimensional periodic lattice of dimensions  = 128 × 128 for 2 × 10 5 simulation steps and considering short influence/long range movement rule.

Figure 6 :
Figure 6: (Color online) spectra of the social structure for the same patterns of Figure 1 row (II) (  = 0.5).In (a) and (b) the existence of a short wavelength structure is observed (in (b) this peak is labeled by an arrow).

Figure 8 :
Figure 8: (Color online) scheme for different interaction energy contributions in the case of the Moore neighborhood in a two-dimensional regular lattice.++ bond is represented in green, −− bond is drawn in orange, +− is drawn in purple, +0 is drawn in sky-blue, −0 is drawn in pink, and finally 00 bond is not drawn.(a) It shows a vertical interface among two populations +1 and −1.The interface energy appears as a consequence of the asymmetry among the individual in the bulk, 8 individuals of the same type, and the individuals at the interface, 5 individuals of one type and 3 from the other.(b) It shows a general interface with both edges and corners.The reader may check the validity of (D.4) and (4).

Figure 9 :
Figure 9: (a) Plots of the rescaled interface edges   /(1 −   ) versus the rescaled number of corners   /(1 −   ) for all the types: +− are symbolized by • and +0 are symbolized by ◼ and ⧫ for −0 case.The colors represent the five different vacancy concentrations: the cyan markers stand for   = 0.01, the red ones stand for   = 0.25, the green ones are for   = 0.5, blue ones stand for   = 0.75, and orange ones stand for   = 0.98.One notices that the correlation is linear:   ≈   .The numerical simulation considers  = 128 2 periodic lattice and includes the 45 different cases and 5 different initial fractions of vacancies.(b) plots the total rescaled number of edges  +− /(1 −   ) versus (,  1 ); (c) the total rescaled number of edges  +− /(1 −   ) versus (,  1 ); and finally (d) plots the total rescaled energy /8(1 −   ) versus (,  1 ).In plots (b), (c), and (d) the surfaces were obtained after an interpolation with a quadratic fit of the form  +  1 +  2 1 +  +  1  +  2 , where a, b, c, d, e, and f are fitting constants.The simulations were developed in a two-dimensional periodic lattice of size  = 128 × 128, considering short influence/long movement range evolution rule and with a running time of 5 × 10 5 steps.

Figure 9 (
d) also presents the relation of the attractor energies as a function of the parameters  and  1 , for the 45 -matrices.A dependence of the rescaled energy as a function of  for various values of  2 is observed.From the data of Figure9(d), it is possible to conclude that the behavior of the corners and the interface length scale are proportional to (1 −   ); that is, / ∼ (1 −   )Φ 1 (,  1 ), where Φ 1 (,  1 ) is a function of  &  1 .Similar relations hold for the number of edges and corners, hence the interface length ℓ  ∼ (1 −   )Φ  (,  1 ).The scaling for the energy  ∼ (1 −   ) follows directly after these relations: because the "bulk" energy becomes  0 = −4( 11  + +  22  − ) = −2(1 −   ) 1 in the special case of  ± = (1 −   )/2, the perimeters ℓ  scale as (1 −   ), hence yielding the total energy scales as  ∼ (1 −   ).