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

Authors: Pablo Medina; Eric Goles; Roberto Zarama; Sergio Rica
Journal: Complexity
 (2017)

---

## Abstract

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.


---

## Body

## 1. 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–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–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 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 q-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. Section2 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.


## 2. The Basic Model

The Sakoda model consists in a regular and periodic lattice withN nodes in a two-dimensional space. Each node i is associated with a discrete variable x
i 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 i is empty (
x
i
=
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 a2
×
2 matrix. We call hereafter the “S-matrix,” which for the sake of space we will symbolize in a vectorized form: S
=
(
s
11
,
s
12
,
s
21
,
s
22
). The entries s
α
β 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 “S-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 individuali:(1)
f
i
x
i
=
∑
k
=
1
N
J
i
k
δ
S
x
i
,
x
k
,where J
i
k
≥
0 denotes a symmetric (J
i
k
=
J
k
i) interaction strength. We avoid self-interactions taking J
i
i
=
0. On the other hand, the sign of the interaction is given by (do not confuse the variable x
i that takes the values x
i
=
0
,
±
1 with the entries of the matrix s
α
β that take the same values s
α
β
=
0
,
±
1):(2)
δ
S
x
i
,
x
k
=
s
α
β
for
x
i
=
α
≠
0
&
x
k
=
β
≠
0
0
o
t
h
e
r
w
i
s
e
.In general, this expression is not symmetric; that is,δ
S
(
x
i
,
x
k
)
≠
δ
S
(
x
k
,
x
i
). We underline that though J
i
k is symmetric, the full interaction J
i
k
δ
S
(
x
i
,
x
k
) is not necessarily symmetric. The coefficients J
i
k 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 i would move towards an empty node j, if f
i
(
x
i
)
<
f
j
(
x
i
). Here f
j
(
x
i
) represents the potential value (1) evaluated at the empty node, j, occupied by the individual x
i. As a general rule, the individual chooses the place j everywhere in the lattice that makes the social expectation function f
j
(
x
i
) the highest. In case of degeneracy, that is, if there are n nodes, j
1
,
j
2
,
…
j
n, such that f
j
1
x
i
=
f
j
2
x
i
=
⋯
=
f
j
n
(
x
i
), 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, ifN
+, N
-, and N
V 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 ϕ
±
=
N
±
/
N and ϕ
V
=
N
V
/
N (naturally ϕ
+
+
ϕ
-
+
ϕ
V
=
1). Hereafter, as in Sakoda original paper [11], we limit our study to the case of ϕ
+
=
ϕ
-
=
(
1
-
ϕ
V
)
/
2.


### 2.1. 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,J
i
k
=
1, if k is in the 8 closest neighbors of the site i, and J
i
k
=
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 possibleS-matrices were run in a two-dimensional periodic lattice of N
=
128
×
128. The evolution usually runs for 2
×
1
0
5 time steps, which happens to be sufficient to reach equilibrium or a steady state. We have also simulated the case of N
=
64
×
64 without significant dependence on the system size.The initial condition consists of a random distribution of states characterized by a fractionϕ
V 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 (ϕ
V
=
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.Figure1 summarizes some social structures resulting from the afore-mentioned algorithm for four characteristic S-matrices and for three different values of the vacancy fraction, ϕ
V. Over all the manuscript, the light gray cells (red cells online) represent individuals at the state +
1 and the dark gray cells (blue cells online) represent -
1 individuals, while the white ones represent the empty spaces, 0.Figure 1
(Color online) snapshots of the few basic social structures: (a)S
=
(
0,1
,
1,0
), (b) S
=
(
1,1
,
-
1
,
-
1
), (c) S
=
(
1,1
,
-
1,1
), and (d) S
=
(
1
,
-
1
,
-
1,1
) for various density of vacancies: (I) ϕ
V
=
0.01, (II) ϕ
V
=
0.5, and (III) ϕ
V
=
0.98. The simulations consider the case of short range of influence and long range of movement and were run in a lattice of N
=
128
×
128 sites and for 2
×
1
0
5 steps.
(a)
(b)
(c)
(d)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 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 ϕ
V
=
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
i
would move towards an empty node
j
, with a probability
p
j
(
x
i
)
=
1
/
Z
e
β
(
f
j
(
x
i
)
-
f
i
(
x
i
)
)
. Here Z
=
1
+
∑
j
‍
e
β
(
f
j
(
x
i
)
-
f
i
(
x
i
)
). 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 self-organize into any spatial structure. Therefore, it is expected that, for a given S-matrix, a critical parameter β
c exists, which controls the existence or not of a social structure.


## 2.1. 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,J
i
k
=
1, if k is in the 8 closest neighbors of the site i, and J
i
k
=
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 possibleS-matrices were run in a two-dimensional periodic lattice of N
=
128
×
128. The evolution usually runs for 2
×
1
0
5 time steps, which happens to be sufficient to reach equilibrium or a steady state. We have also simulated the case of N
=
64
×
64 without significant dependence on the system size.The initial condition consists of a random distribution of states characterized by a fractionϕ
V 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 (ϕ
V
=
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.Figure1 summarizes some social structures resulting from the afore-mentioned algorithm for four characteristic S-matrices and for three different values of the vacancy fraction, ϕ
V. Over all the manuscript, the light gray cells (red cells online) represent individuals at the state +
1 and the dark gray cells (blue cells online) represent -
1 individuals, while the white ones represent the empty spaces, 0.Figure 1
(Color online) snapshots of the few basic social structures: (a)S
=
(
0,1
,
1,0
), (b) S
=
(
1,1
,
-
1
,
-
1
), (c) S
=
(
1,1
,
-
1,1
), and (d) S
=
(
1
,
-
1
,
-
1,1
) for various density of vacancies: (I) ϕ
V
=
0.01, (II) ϕ
V
=
0.5, and (III) ϕ
V
=
0.98. The simulations consider the case of short range of influence and long range of movement and were run in a lattice of N
=
128
×
128 sites and for 2
×
1
0
5 steps.
(a)
(b)
(c)
(d)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 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 ϕ
V
=
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
i
would move towards an empty node
j
, with a probability
p
j
(
x
i
)
=
1
/
Z
e
β
(
f
j
(
x
i
)
-
f
i
(
x
i
)
)
. Here Z
=
1
+
∑
j
‍
e
β
(
f
j
(
x
i
)
-
f
i
(
x
i
)
). 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 self-organize into any spatial structure. Therefore, it is expected that, for a given S-matrix, a critical parameter β
c exists, which controls the existence or not of a social structure.


## 3. Potts Energy for the Sakoda Model

### 3.1. 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 theS-matrix is symmetric, then the Potts-like energy function [24, 25], of a given configuration of the network x
=
{
x
1
,
x
2
,
…
,
x
i
,
…
,
x
k
,
…
},(3)
E
x
=
-
1
2
∑
i
,
k
J
i
k
δ
S
x
i
,
x
kis a nonincreasing function after a movement resulting from the Sakoda algorithm; that is, E
a
f
t
e
r
≤
E
b
e
f
o
r
e. In other words, E
[
x
] is a nonstrict Lyapunov function.Proof.
Suppose thatS is a symmetric matrix. Thus, δ
S
(
x
j
,
x
k
)
=
δ
S
(
x
k
,
x
j
) for two arbitrary sites j and k in the lattice. Consider the exchange of the sites l (occupied) and q (empty), that is, x
l
=
±
1 and x
q
=
0. Hence, the initial x and the subsequent configuration x
′ after a swap between the q and l sites are x
=
{
…
,
x
l
,
…
,
x
q
=
0
,
…
} and x
′
=
{
…
,
x
q
=
0
,
…
,
x
l
,
…
}
. If the S-matrix is symmetric, energy (3) difference between the initial and the actual configuration becomes(4)
Δ
E
=
E
x
′
-
E
x
=
∑
j
≠
q
,
l
J
l
j
-
J
q
j
δ
S
x
l
,
x
j
.Here we have used that the self-interaction, J
l
l
=
0, vanishes and δ
S
(
x
q
,
x
j
)
=
0 for all j. This occurs because the δ
S-function vanishes exactly if it is evaluated at a vacancy state x
q
=
0 (for details see Appendix B). The intrinsic mechanism of Sakoda’s rule imposes that the social expectation (1), f
l, of the individual x
l, is smaller than the social expectation at the site q, that is, f
q
(
x
l
)
>
f
l
(
x
l
):(5)
f
q
x
l
=
∑
j
J
q
j
δ
S
x
l
,
x
j
>
f
l
x
l
=
∑
j
J
l
j
δ
S
x
l
,
x
j
.Rewriting the r.h.s. of (4) in terms of the expressions f
l and f
q of (1) and (5), respectively, we obtain Δ
E
=
f
l
(
x
l
)
-
f
q
(
x
l
)
≤
0
. We emphasize that the equality holds in the case of f
q
(
x
l
)
=
f
l
(
x
l
); in this case no exchange is carried out. In conclusion, if S is a symmetric matrix, then Δ
E
≤
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 theS-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 F
=
E
-
μ
+
N
+
-
μ
-
N
-
-
μ
V
N
V, where μ
± and μ
V 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).In AppendixC, we show that the energy (3) corresponds, in the present case, exactly to the spin 1 Ising model [20]. Moreover, if the S-matrix is S
=
(
1
,
-
1
,
-
1,1
), one readily gets E
[
x
]
=
-
1
/
2
∑
i
,
k
J
i
k
x
i
x
k, which is exactly the expression of the Ising energy associated with Schelling’s model [21].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(6)
E
=
-
4
s
11
N
+
+
s
22
N
-
+
3
2
χ
l
+
-
+
s
11
l
+
0
+
s
22
l
-
0
,
χ
=
s
11
+
s
22
-
s
12
+
s
21
,where l
α
β
=
p
α
β
-
1
/
3
e
α
β is an effective interface length which is directly related to the number of edges, p
α
β, and corners, e
α
β, among the interface created by the states α and β (see Appendix D for a derivation of (6)). Finally, the term -
4
(
s
11
N
+
+
s
22
N
-
) represents the bulk energy contribution. However, because N
±
=
N
ϕ
± 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.


### 3.2. Energy Behavior in the Sakoda Model

To illustrate the energy principle, Figure2(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).Figure 2
(Color online) temporal evolution of the energy (i) and the effective perimetersl
α
β (ii) for the same cases as in Figure 1 row (II), that is, for ϕ
V
=
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 S-matrix; hence its energy may decrease or increase indistinctly. Some increments of energy are labeled with ∗. (ii) The perimeters l
α
β 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
×
1
0
6 time steps.Quantitatively, the energy expansion (6) allows us to distinguish between the segregation and aggregation patterns. Indeed, the energy contributions E
+
-
=
3
/
2
χ
l
+
-, E
+
0
=
3
/
2
s
11
l
+
0, and E
-
0
=
3
/
2
s
22
l
-
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 χ, s
11, and s
22 would be a kind of “social tension” among the interfaces. From a physical point of view, an interface isnonmiscible if the respective social tension, namely, χ, s
11, or s
22, is positive. Hence, because the system minimizes the energy, it minimizes the corresponding perimeter; otherwise, the interface would bemiscible. 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 Figure1(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 Figure1 may be understood according to (6). In Figure 1(a), that is, S
=
(
0,1
,
1,0
) and χ
=
-
2, the perimeter energy reads -
3
l
+
-. In this case, the energy decreases if the interface perimeter between -
1 and +
1 populations increases (see l
+
- blue curve in Figure 2(b)), thus creating a labyrinth pattern in analogy with patterns of antiferromagnetic materials.In the case of Figure1(b), that is, S
=
(
1,1
,
-
1
,
-
1
) and χ
=
0, the energy term is simplified to 3
/
2
(
l
+
0
-
l
-
0
), being this a minima if, simultaneously, l
+
0 is maximized (see the upper red curve in Figure 2(b)) and l
-
0 is minimized (see the l
+
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), S
=
(
1,1
,
-
1,1
), that is, χ
=
-
2, and Figure 1(d), S
=
(
1
,
-
1
,
-
1,1
) (χ
=
4), produce segregation patterns. For these cases, the perimeter energy terms are similar: 1
/
2
(
2
l
+
-
+
l
+
0
+
l
-
0
) and 3
/
2
(
4
l
+
-
+
l
+
0
+
l
-
0
) for Figures 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 l
+
- (d)), the interface perimeter, l
+
-, 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.


## 3.1. 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 theS-matrix is symmetric, then the Potts-like energy function [24, 25], of a given configuration of the network x
=
{
x
1
,
x
2
,
…
,
x
i
,
…
,
x
k
,
…
},(3)
E
x
=
-
1
2
∑
i
,
k
J
i
k
δ
S
x
i
,
x
kis a nonincreasing function after a movement resulting from the Sakoda algorithm; that is, E
a
f
t
e
r
≤
E
b
e
f
o
r
e. In other words, E
[
x
] is a nonstrict Lyapunov function.Proof.
Suppose thatS is a symmetric matrix. Thus, δ
S
(
x
j
,
x
k
)
=
δ
S
(
x
k
,
x
j
) for two arbitrary sites j and k in the lattice. Consider the exchange of the sites l (occupied) and q (empty), that is, x
l
=
±
1 and x
q
=
0. Hence, the initial x and the subsequent configuration x
′ after a swap between the q and l sites are x
=
{
…
,
x
l
,
…
,
x
q
=
0
,
…
} and x
′
=
{
…
,
x
q
=
0
,
…
,
x
l
,
…
}
. If the S-matrix is symmetric, energy (3) difference between the initial and the actual configuration becomes(4)
Δ
E
=
E
x
′
-
E
x
=
∑
j
≠
q
,
l
J
l
j
-
J
q
j
δ
S
x
l
,
x
j
.Here we have used that the self-interaction, J
l
l
=
0, vanishes and δ
S
(
x
q
,
x
j
)
=
0 for all j. This occurs because the δ
S-function vanishes exactly if it is evaluated at a vacancy state x
q
=
0 (for details see Appendix B). The intrinsic mechanism of Sakoda’s rule imposes that the social expectation (1), f
l, of the individual x
l, is smaller than the social expectation at the site q, that is, f
q
(
x
l
)
>
f
l
(
x
l
):(5)
f
q
x
l
=
∑
j
J
q
j
δ
S
x
l
,
x
j
>
f
l
x
l
=
∑
j
J
l
j
δ
S
x
l
,
x
j
.Rewriting the r.h.s. of (4) in terms of the expressions f
l and f
q of (1) and (5), respectively, we obtain Δ
E
=
f
l
(
x
l
)
-
f
q
(
x
l
)
≤
0
. We emphasize that the equality holds in the case of f
q
(
x
l
)
=
f
l
(
x
l
); in this case no exchange is carried out. In conclusion, if S is a symmetric matrix, then Δ
E
≤
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 theS-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 F
=
E
-
μ
+
N
+
-
μ
-
N
-
-
μ
V
N
V, where μ
± and μ
V 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).In AppendixC, we show that the energy (3) corresponds, in the present case, exactly to the spin 1 Ising model [20]. Moreover, if the S-matrix is S
=
(
1
,
-
1
,
-
1,1
), one readily gets E
[
x
]
=
-
1
/
2
∑
i
,
k
J
i
k
x
i
x
k, which is exactly the expression of the Ising energy associated with Schelling’s model [21].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(6)
E
=
-
4
s
11
N
+
+
s
22
N
-
+
3
2
χ
l
+
-
+
s
11
l
+
0
+
s
22
l
-
0
,
χ
=
s
11
+
s
22
-
s
12
+
s
21
,where l
α
β
=
p
α
β
-
1
/
3
e
α
β is an effective interface length which is directly related to the number of edges, p
α
β, and corners, e
α
β, among the interface created by the states α and β (see Appendix D for a derivation of (6)). Finally, the term -
4
(
s
11
N
+
+
s
22
N
-
) represents the bulk energy contribution. However, because N
±
=
N
ϕ
± 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.


## 3.2. Energy Behavior in the Sakoda Model

To illustrate the energy principle, Figure2(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).Figure 2
(Color online) temporal evolution of the energy (i) and the effective perimetersl
α
β (ii) for the same cases as in Figure 1 row (II), that is, for ϕ
V
=
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 S-matrix; hence its energy may decrease or increase indistinctly. Some increments of energy are labeled with ∗. (ii) The perimeters l
α
β 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
×
1
0
6 time steps.Quantitatively, the energy expansion (6) allows us to distinguish between the segregation and aggregation patterns. Indeed, the energy contributions E
+
-
=
3
/
2
χ
l
+
-, E
+
0
=
3
/
2
s
11
l
+
0, and E
-
0
=
3
/
2
s
22
l
-
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 χ, s
11, and s
22 would be a kind of “social tension” among the interfaces. From a physical point of view, an interface isnonmiscible if the respective social tension, namely, χ, s
11, or s
22, is positive. Hence, because the system minimizes the energy, it minimizes the corresponding perimeter; otherwise, the interface would bemiscible. 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 Figure1(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 Figure1 may be understood according to (6). In Figure 1(a), that is, S
=
(
0,1
,
1,0
) and χ
=
-
2, the perimeter energy reads -
3
l
+
-. In this case, the energy decreases if the interface perimeter between -
1 and +
1 populations increases (see l
+
- blue curve in Figure 2(b)), thus creating a labyrinth pattern in analogy with patterns of antiferromagnetic materials.In the case of Figure1(b), that is, S
=
(
1,1
,
-
1
,
-
1
) and χ
=
0, the energy term is simplified to 3
/
2
(
l
+
0
-
l
-
0
), being this a minima if, simultaneously, l
+
0 is maximized (see the upper red curve in Figure 2(b)) and l
-
0 is minimized (see the l
+
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), S
=
(
1,1
,
-
1,1
), that is, χ
=
-
2, and Figure 1(d), S
=
(
1
,
-
1
,
-
1,1
) (χ
=
4), produce segregation patterns. For these cases, the perimeter energy terms are similar: 1
/
2
(
2
l
+
-
+
l
+
0
+
l
-
0
) and 3
/
2
(
4
l
+
-
+
l
+
0
+
l
-
0
) for Figures 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 l
+
- (d)), the interface perimeter, l
+
-, 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.


## 4. The Basic Social Structures

Under the same conditions of Figure1, we have run numerical simulations for all 45 possible S-matrices. The social structures may be classified in terms of the values of the S-matrix based on three different aspects: (1) the particular social structure, that is, the pattern, and (2) the pattern power spectrum which provides information on the spatial correlations (see Section 5.3). The best classification scheme appears to be in terms of the parameters ζ
1
=
s
11
+
s
22 and ζ
2
=
s
12
+
s
21. The first parameter gives the information of the attitudes towards their own group, while the second provides information about the attitudes of individuals towards the opposite group. Considering ζ
1 and ζ
2, it makes classifying 25 different cases possible as shown in Figure 3 and in Table 1.Table 1
Summary information on the patterns of Figure3. 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.
ζ
1
−2
−1
0
1
2
ζ
2
2
S
=
(
-
1,1
,
1
,
-
1
)
S
=
(
0,1
,
1
,
-
1
)
S
=
(
1,1
,
1
,
-
1
)
S
=
(
1,1
,
1,0
)
S
=
(
1,1
,
1,1
)
μ = none
μ
=
-
0.5
μ
=
-
1.4
μ
=
-
0.3
μ = none
1
S
=
(
-
1,1
,
0
,
-
1
)
S
=
(
0,1
,
0
,
-
1
)
S
=
(
1,1
,
0
,
-
1
)
S
=
(
1,1
,
0,0
)
S
=
(
1,1
,
0,1
)
μ
=
n
o
n
e
μ
=
-
0.4
μ
=
-
2.5
μ
=
-
1.2
μ
=
-
3
0
S
=
(
-
1,1
,
-
1
,
-
1
)
S
=
(
0,1
,
-
1
,
-
1
)
S
=
(
1,1
,
-
1
,
-
1
)
S
=
(
1,1
,
-
1,0
)
S
=
(
1,0
,
0,1
)
μ = none
μ = none
μ
=
-
3
μ
=
-
1.6
μ
=
-
3.5
−1
S
=
(
-
1,0
,
-
1
,
-
1
)
S
=
(
0,0
,
-
1
,
-
1
)
S
=
(
1,0
,
-
1
,
-
1
)
S
=
(
1,0
,
-
1,0
)
S
=
(
1,0
,
-
1,1
)
μ = none
μ = none
μ
=
-
1.4
μ
=
-
2.5
μ
=
-
4
−2
S
=
(
-
1
,
-
1
,
-
1
,
-
1
)
S
=
(
0
,
-
1
,
-
1
,
-
1
)
S
=
(
0
,
-
1
,
-
1,0)
S
=
(
1
,
-
1
,
-
1,0
)
S
=
(
1
,
-
1
,
-
1,1
)
μ = none
μ
=
-
2.2
μ
=
-
3
μ
=
-
3
μ
=
-
4Figure 3
(Color online) classification of the attractors for differentS-matrices according to the ζ
1
=
s
11
+
s
22 and ζ
2
=
s
12
+
s
21 plane. The precise S-matrices are given in Table 1. The 25 numerical simulations were done in a 128
×
128 lattice with random initial configurations and ϕ
V
=
1
/
2. As expected from the energy expansion (6), the diagonal s
11
+
s
22
=
s
12
+
s
21 splits roughly the segregation-like and the aggregation-like patterns: +
1/-
1 segregation behavior is observed below this line, while a heterogeneous attraction behavior is mostly observed above this line. We highlight the following social structures:segregation (S1),crossroads (S2),mutual suspicion (S3),social workers (S5),husband and wife (S6),boys and girls (S7),inclusion (S8),individualism (S9), andexclusion (S10), which are representative of the collective behavior.The phase diagram of Figure3 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, andhusband and wife. Based on the symmetries of the diagram, we suggest in the present work three new social behaviors:Inclusion,Repulsion, andExclusion, 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χ
=
s
11
+
s
22
-
(
s
12
+
s
21
) splits roughly the segregation-like and the aggregation-like pattern.The description of the social behaviors presented in Figure3 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)).
(S2)
Crossroads manifests also segregation; however in this case, the individuals of other groups do not interact. The S-matrix reads S
=
(
1,0
,
0,1
)
.
(S3)
Mutual suspicionis a characteristic typified by neutrality among members of the same group and repulsion towards the members of the other group. In this case the S-matrix reads S
=
(
0
,
-
1
,
-
1,0
).
(S4)
Social climbers represents the interaction between a lower class and an upper class this is characterized by a S-matrix S
=
(
1
,
-
1,1
,
-
1
).
(S5)
Social workers is represented in Figure 1(b). In this case, one group feels attraction by everybody while the other group feels repulsion for everybody. Its S-matrix reads S
=
(
1,1
,
-
1
,
-
1
).
(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 S
=
(
1,1
,
1
,
-
1
).
(S7)
Boys and girls is exactly the opposite of segregation, S
=
(
-
1,1
,
1
,
-
1
). This is an example of antiferromagnetic physics. In this case (χ
=
-
4), the perimeter is maximized [22].
(S8)
Inclusion manifests the social situations where all individuals want to be always together. In this case S
=
(
1,1
,
1,1
).
(S9)
Individualism is a case where the individuals want to be as far as possible from others. Here S
=
(
-
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 l
+
-, 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 S-matrix is S
=
(
-
1,1
,
-
1
,
-
1
).


## 5. Quantitative Description of the Social Structures

### 5.1. Phase Diagram in theχ-E 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
-
ϕ
V
)
/
2, the energy expression (6) reads (N
±
=
ϕ
±
N) (7)
E
=
-
2
N
1
-
ϕ
V
s
11
+
s
22
+
3
2
χ
l
+
-
+
s
11
l
+
0
+
s
22
l
-
0
. From the numerical simulations, we conclude that the behavior of the corners and the interface length scale proportional also to (
1
-
ϕ
V
), that is, the total energy scales, as E
/
N
~
(
1
-
ϕ
V
)
Φ
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.Figure4 shows a set of 40 of the 45 attractors (the missing 5 social structures are overlapped) classified in the χ-E plane for ϕ
V
=
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 S
=
(
-
1
,
-
1
,
-
1
,
-
1
) that correspond to χ
=
0 and ζ
1
=
s
11
+
s
22
=
-
2 and an absolute minimum S
=
(
1,1
,
1,1
) (χ
=
0, ζ
1
=
2). The social structures are spanned all over the χ-E plane below this envelope (
i
) and above the line (
i
i
) for χ
<
0 and the horizontal line E
/
4
(
1
-
ϕ
V
)
N
=
-
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
(Color online) phase diagram of 40 attractors inχ
-
E plane (normalized energy per occupied site, E
/
4
(
1
-
ϕ
V
)
N) for ϕ
V
=
0.5 (notice that the attractors (
0
,
-
1,1
,
-
1
), (
1,0
,
0
,
-
1
), (
0,1
,
-
1,0
), (
0
,
-
1,0
,
-
1
), and (
1,1
,
-
1,0
) are overlapped by other attractors in the figure). These simulations correspond to a N
=
128
×
128 two-dimensional lattice running up to 5
×
1
0
5 simulation steps in the case of short influence/long range movement evolution rule. Attractors with social interpretations are labeled (from (S1) to (S11)). Notice the patternsocial climber (S4) missing in Figure 3. The missing five attractors are overlapped by others because of an approximate degeneracy of the energy for the same parameter χ. These are (
0
,
-
1,1
,
-
1
), (
1,0
,
0
,
-
1
), (
0,1
,
-
1,0
), (
0
,
-
1,0
,
-
1
), and (
1,1
,
-
1,0
). The curves (
i
) and (
i
i
) are sketches of the upper and lower bounds, respectively (see text).Figure4 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 (
E
,
χ
) 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 asNeutrals. The social structures belong toNeutrals class where the individuals of different groups do not express neither rejection nor attraction among others. This behavior is represented by the interaction S
=
(
0,0
,
0,0
). Similarly, we identify the 9th group of social behavior (G9) named bysocial neutrality.


### 5.2. The Role of the Vacancy Fraction on the Social Structures

As it can be noticed in Figure1, the final patterns depend on the empty spaces in the lattice, ϕ
V. In particular, we observe that the number of representative social structures diminishes as both ϕ
V
→
0 and also in the limit ϕ
V
→
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 ϕ
V on the final patterns.Figure5 summarizes the phase diagram of some attractors in the χ
-
ϕ
V space. The reader may note that different vacancy concentrations, ϕ
V, 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, S
=
(
1
,
-
1
,
-
1,1
), one observes that if we increase the number of total vacancies we recover the phenomenology already found in [20].Figure 5
(Color online) examples of different attractors classified byχ parameter for different values of ϕ
V as indicated in the plot. Numerical simulations run in a two-dimensional periodic lattice of dimensions N
=
128
×
128 for 2
×
1
0
5 simulation steps and considering short influence/long range movement rule.


### 5.3. Existence of Long Range Order on Social Structures

Generally speaking, it is observed as a common rule that all fractionsϕ
V present phases of segregation (ferromagnetic phase), heterogeneous attraction (antiferromagnetic phase), and a neutral or disordered phase (paramagnetic phase). More deeply, in a close analogy with phase transitions, we observe that the segregation-like patterns do manifest long range order (more precisely quasi-long range order), but the cases of heterogeneous attraction and disordered phases do not.The quasi-long range order is readily quantified by the spatial correlation function which readsC
(
r
)
=
1
/
N
∑
l
x
l
x
l
+
r. It may be computed directly in terms of the Fourier transform of the spatial pattern: (8)
x
^
q
=
1
N
∑
l
x
l
e
2
π
i
q
l
/
L
.Here x
l is the spatial pattern, L
=
N is the linear size of a square lattice, and q is an integer ranging within -
L
/
2
+
1
≤
q
≤
L
/
2. A straight-forward calculation relates the correlation function with the spectrum. In the limit L
,
N
→
∞, we have (9)
x
k
=
∑
q
x
^
q
e
-
2
π
i
q
k
/
L
,thus yielding(10)
C
r
=
∑
q
,
p
x
^
q
x
^
p
e
2
π
i
q
r
/
L
1
N
∑
l
e
2
π
i
q
+
p
l
L
⟶
∑
q
x
^
q
2
e
2
π
i
q
r
/
L
. 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 q): I
q
=
x
^
q
2.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 well-defined domains, one expects thatI
q
~
q
-
μ, 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 ϕ
V
=
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 q) 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.Figure 6
(Color online) spectra of the social structure for the same patterns of Figure1 row (II) (ϕ
V
=
0.5). In (a) and (b) the existence of a short wavelength structure is observed (in (b) this peak is labeled by an arrow).
(a)
(b)
(c)
(d)


## 5.1. Phase Diagram in theχ-E 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
-
ϕ
V
)
/
2, the energy expression (6) reads (N
±
=
ϕ
±
N) (7)
E
=
-
2
N
1
-
ϕ
V
s
11
+
s
22
+
3
2
χ
l
+
-
+
s
11
l
+
0
+
s
22
l
-
0
. From the numerical simulations, we conclude that the behavior of the corners and the interface length scale proportional also to (
1
-
ϕ
V
), that is, the total energy scales, as E
/
N
~
(
1
-
ϕ
V
)
Φ
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.Figure4 shows a set of 40 of the 45 attractors (the missing 5 social structures are overlapped) classified in the χ-E plane for ϕ
V
=
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 S
=
(
-
1
,
-
1
,
-
1
,
-
1
) that correspond to χ
=
0 and ζ
1
=
s
11
+
s
22
=
-
2 and an absolute minimum S
=
(
1,1
,
1,1
) (χ
=
0, ζ
1
=
2). The social structures are spanned all over the χ-E plane below this envelope (
i
) and above the line (
i
i
) for χ
<
0 and the horizontal line E
/
4
(
1
-
ϕ
V
)
N
=
-
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
(Color online) phase diagram of 40 attractors inχ
-
E plane (normalized energy per occupied site, E
/
4
(
1
-
ϕ
V
)
N) for ϕ
V
=
0.5 (notice that the attractors (
0
,
-
1,1
,
-
1
), (
1,0
,
0
,
-
1
), (
0,1
,
-
1,0
), (
0
,
-
1,0
,
-
1
), and (
1,1
,
-
1,0
) are overlapped by other attractors in the figure). These simulations correspond to a N
=
128
×
128 two-dimensional lattice running up to 5
×
1
0
5 simulation steps in the case of short influence/long range movement evolution rule. Attractors with social interpretations are labeled (from (S1) to (S11)). Notice the patternsocial climber (S4) missing in Figure 3. The missing five attractors are overlapped by others because of an approximate degeneracy of the energy for the same parameter χ. These are (
0
,
-
1,1
,
-
1
), (
1,0
,
0
,
-
1
), (
0,1
,
-
1,0
), (
0
,
-
1,0
,
-
1
), and (
1,1
,
-
1,0
). The curves (
i
) and (
i
i
) are sketches of the upper and lower bounds, respectively (see text).Figure4 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 (
E
,
χ
) 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 asNeutrals. The social structures belong toNeutrals class where the individuals of different groups do not express neither rejection nor attraction among others. This behavior is represented by the interaction S
=
(
0,0
,
0,0
). Similarly, we identify the 9th group of social behavior (G9) named bysocial neutrality.


## 5.2. The Role of the Vacancy Fraction on the Social Structures

As it can be noticed in Figure1, the final patterns depend on the empty spaces in the lattice, ϕ
V. In particular, we observe that the number of representative social structures diminishes as both ϕ
V
→
0 and also in the limit ϕ
V
→
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 ϕ
V on the final patterns.Figure5 summarizes the phase diagram of some attractors in the χ
-
ϕ
V space. The reader may note that different vacancy concentrations, ϕ
V, 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, S
=
(
1
,
-
1
,
-
1,1
), one observes that if we increase the number of total vacancies we recover the phenomenology already found in [20].Figure 5
(Color online) examples of different attractors classified byχ parameter for different values of ϕ
V as indicated in the plot. Numerical simulations run in a two-dimensional periodic lattice of dimensions N
=
128
×
128 for 2
×
1
0
5 simulation steps and considering short influence/long range movement rule.


## 5.3. Existence of Long Range Order on Social Structures

Generally speaking, it is observed as a common rule that all fractionsϕ
V present phases of segregation (ferromagnetic phase), heterogeneous attraction (antiferromagnetic phase), and a neutral or disordered phase (paramagnetic phase). More deeply, in a close analogy with phase transitions, we observe that the segregation-like patterns do manifest long range order (more precisely quasi-long range order), but the cases of heterogeneous attraction and disordered phases do not.The quasi-long range order is readily quantified by the spatial correlation function which readsC
(
r
)
=
1
/
N
∑
l
x
l
x
l
+
r. It may be computed directly in terms of the Fourier transform of the spatial pattern: (8)
x
^
q
=
1
N
∑
l
x
l
e
2
π
i
q
l
/
L
.Here x
l is the spatial pattern, L
=
N is the linear size of a square lattice, and q is an integer ranging within -
L
/
2
+
1
≤
q
≤
L
/
2. A straight-forward calculation relates the correlation function with the spectrum. In the limit L
,
N
→
∞, we have (9)
x
k
=
∑
q
x
^
q
e
-
2
π
i
q
k
/
L
,thus yielding(10)
C
r
=
∑
q
,
p
x
^
q
x
^
p
e
2
π
i
q
r
/
L
1
N
∑
l
e
2
π
i
q
+
p
l
L
⟶
∑
q
x
^
q
2
e
2
π
i
q
r
/
L
. 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 q): I
q
=
x
^
q
2.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 well-defined domains, one expects thatI
q
~
q
-
μ, 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 ϕ
V
=
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 q) 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.Figure 6
(Color online) spectra of the social structure for the same patterns of Figure1 row (II) (ϕ
V
=
0.5). In (a) and (b) the existence of a short wavelength structure is observed (in (b) this peak is labeled by an arrow).
(a)
(b)
(c)
(d)


## 6. 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 parameterss
11 and s
22 and χ
=
(
s
11
+
s
22
)
-
(
s
12
+
s
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 Figures3 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 S-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.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 theS-matrix for the interactions among states and vacancies. In this case, the S-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 q-states of the variable x
i, the S-matrix becomes q
×
q interaction matrix and the total possible number of different cases to be considered would become huge: 3
q
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.Figure7 summarizes the attractors for the same S-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 pattern depends on the range of interaction and/or the scale of the movement. In those attractors that exhibit patterns of aggregation or segregation, long movements produce larger clusters of individuals than short movements. The range of interaction also plays a role on the scale of the patterns in the social structures. The largest pattern arises for the case of long range movement long range interaction (Figure 7(III)); the next intermediate size would be for the case of short range interaction long range movement (Figure 1(II)). For the cases of short range movement, the pattern’s scale is usually smaller. Figure 7(I) shows that the case of short range interaction-short range movement displays the smallest pattern’s size. Finally, Figure 7(II) corresponds to the case of long range interaction-short range movement with a slightly larger scale pattern. 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.Figure 7
Snapshots for different interaction and movement (rest) rules for the same cases as in Figure1. (a) S
=
(
0,1
,
1,0
), (b) S
=
(
1,1
,
-
1
,
-
1
), (c) S
=
(
1,1
,
-
1,1
), and (d) S
=
(
1
,
-
1
,
-
1,1
). The first line (I) corresponds to the case of short range interaction and short range movement. The second row (II) corresponds to long range interaction/short range movement, and, finally the third line (III) corresponds to the case of long range interaction and long range movement. Simulations were run in a N
=
128
×
128 lattice with ϕ
V
=
0.5 and for 2
×
1
0
5 steps for (I) and (II) and 5
×
1
0
4 steps for (III).
(a)
(b)
(c)
(d)



---

*Source: 3548591-2017-01-23.xml*