# Socioeconomic agents as active matter in nonequilibrium Sakoda-Schelling
models

Ruben Zakine<EMAIL_ADDRESS>Chair of Econophysics and Complex
Systems, École polytechnique, 91128 Palaiseau Cedex, France  LadHyX, CNRS,
École polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France
Jérôme Garnier-Brun  Chair of Econophysics and Complex Systems, École
polytechnique, 91128 Palaiseau Cedex, France  LadHyX, CNRS, École
polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France
Antoine-Cyrus Becharat  Chair of Econophysics and Complex Systems, École
polytechnique, 91128 Palaiseau Cedex, France  LadHyX, CNRS, École
polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France
Michael Benzaquen  Chair of Econophysics and Complex Systems, École
polytechnique, 91128 Palaiseau Cedex, France  LadHyX, CNRS, École
polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France
Capital Fund Management, 23 Rue de l’Université, 75007 Paris, France

(September 3, 2024)

###### Abstract

How robust are socioeconomic agent-based models with respect to the details of
the agents’ decision rule? We tackle this question by considering an
occupation model in the spirit of the Sakoda-Schelling model, historically
introduced to shed light on segregation dynamics among human groups. For a
large class of utility functions and decision rules, we pinpoint the
nonequilibrium nature of the agent dynamics, while recovering the equilibrium-
like phase separation phenomenology. Within the mean-field approximation we
show how the model can be mapped, to some extent, onto an active matter field
description. Finally, we consider non-reciprocal interactions between two
populations, and show how they can lead to non-steady macroscopic behavior. We
believe our approach provides a unifying framework to further study geography-
dependent agent-based models, notably paving the way for joint consideration
of population and price dynamics within a field theoretic approach.

##  I Introduction

Will a collective system in which individuals share the same common goal ever
reach an optimal state? This nontrivial question is at the very core of strong
debates among economists, notably because the notion of “optimal state” is
intrinsically political and most often ill-defined. Despite the common idea
that a system made of agents individualistically improving their outcome will
spontaneously converge by the action of the “invisible hand” to an optimal
collective state, simple models have been shown to contradict this belief [1,
2, 3]. A well documented example of such system is the celebrated Schelling
model [4]. The latter can be considered to be a variant of the model
previously111To be perfectly precise, the first mention of Sakoda’s model can
be traced back to his unpublished PhD thesis completed in 1946, while
Schelling’s work can be found in a 1969 working paper [5]. In any case, there
is no reason to think either author took inspiration from the other, the
objective of the papers being clearly quite different. introduced by Sakoda
[6], and will thus be referred henceforth as the Sakoda-Schelling model. To
understand some aspects of urban segregation in post-WWII American cities, and
more widely of urban and social dynamics, both authors proposed simple lattice
models of idealized cities. Each site, representing an accommodation, can be
empty or occupied by an agent belonging to one of two sub-populations in the
system. Interestingly, Schelling observed that when introducing a slight
preference for agents to be surrounded by neighbors of their own group, the
system evolves towards configurations with completely segregated regions.
While in fact not very well suited to explain urban segregation, which is
intimately related to past and present public policies rather than self-
organization [7, 8], the model illustrates how the micromotives of the agents
may lead to unanticipated macrobehavior [9].

Along the years, the Sakoda-Schelling model has attracted further attention of
statistical physicists [10, 11, 12, 13], due to its simple microscopic rules,
its paradoxical macroscopic consequences and its unconventional non-local
particle moves. To the usual end of bridging the gap from _micro to macro_ ,
mappings onto equilibrium systems were suggested [14], but with limited
analytical results. To gain a more in-depth understanding of the mechanism
through which individual choices may lead to sub-optimal collective outcomes,
Grauwin et al. introduced a modified version of the Schelling model with a
single type of agent occupying a lattice divided in pre-defined neighborhoods,
or blocks [15]. In this occupation model, the agents now base their decisions
on the neighborhood density, which is identical for all the agents in a given
block. This fixed neighborhood structure then allows to describe analytically
the steady state as the minimizer of a free energy, and to recover a
nontrivial phase with suboptimal jam-packed neighborhoods. Subsequent works
have then explored variations of these different models focusing on the effect
of altruistic agents [16], dynamics close to criticality [17, 18, 19] or habit
formation [20].

Even in the seemingly simpler occupation problem of Grauwin et al. [15],
several questions persist, both from the socioeconomic and statistical physics
perspectives. In particular, the role of the specific decision rule and the
precise nature of neighborhoods on the phenomenology of the model remain
unclear. Indeed, to allow for the standard techniques of statistical mechanics
to be applicable, the choice of the neighborhoods and the dynamics is very
constrained, see [21]. As will be discussed in detail, most non-trivial
decision rules lead the system out of thermodynamic equilibrium, requiring
calculations that are not always readily tractable. As it is extremely
difficult to empirically determine how economic agents actually make
decisions, the physics-inspired theoretical analysis of toy models has a
significant part to play, in particular to determine the robustness of
qualitative findings to specific modeling choices. Besides, as argued in [21]
and by some of us in [22], the intrinsically individualistic nature of agent-
specific moves in socioeconomic models means that the description of
collective behaviors as the minimization of some global energy is often not
possible. Understanding simple out-of-equilibrium dynamics as those that arise
from the decision rules presented here is therefore also necessary from the
methodological point of view.

The purpose of this paper is to assess, within a general Sakoda-Schelling like
occupation model, whether and how the sub-optimal concentration of agents in
overly dense regions still occurs out of equilibrium. Most importantly, we
relax the assumption of taking a specific decision rule, and no longer require
pre-defined block neighborhoods as in [15]. The resulting heterogeneity of
interactions in our model then requires the use of out-of-equilibrium
statistical mechanics techniques, the progress of which in the last decade can
be credited to active matter theory. Overall, we find that the phenomenology
of the model is largely unaffected by its nonequilibrium nature, suggesting
that the tendency of agents to aggregate sub-optimally is robust to large
classes of decision rules. This being said, our analysis highlights
interesting theoretical subtleties, notably related to the non-monotonicity of
the utility functions considered, that may, in turn, contribute to the
understanding of other complex physical systems.

The paper is organized as follows. In Sec. II we introduce a Schelling-like
occupation model, in which we keep the utility function and decision rule as
general as possible to allow for nonequilibrium dynamics. We the perform a
numerical analysis of the model. In Sec. III we present a mean-field
description of the dynamics, and determine the region in parameter space where
condensation necessarily occurs. In Sec. IV we show how the dynamics can be
mapped on the Active Model B [23], which is considered to be the natural
nonequilibrium extension of the Cahn-Hilliard field relaxation [24]. This
mapping notably allows to compute the phase densities of the concentrated
states. In Sec. V we propose some relevant generalizations of the model,
namely with two different populations and a housing market. Finally, in Sec.
VI we discuss the implications of our study and conclude.

##  II A Sakoda-Schelling occupation model

###  II.1 Setup

Consider a city structured as a two-dimensional rectangular lattice composed
of $M=L_{x}\times L_{y}$ sites (or houses). Each site can be occupied by at
most one of the $N(\leq M)$ agents living in this city. On each site of
coordinate $\bm{r}=(i,j)$, the occupation field $n$ takes the value
$n(\bm{r})=1$ if the site is occupied, $n(\bm{r})=0$ if it is vacant. It is
assumed that each agent $k$ wants to maximize their own utility $u_{k}$, which
depends on the local density of agents around them. Typically, it is natural
to think that people like to gather in relatively dense areas to benefit from
the city life, but not too dense as crowding might degrade the quality of
life. Agents estimate the local density by averaging the occupation field with
a probability-density-function kernel $G_{\sigma}$, where $\sigma$ stands for
the interaction range. The kernel is assumed to be isotropic and identical for
all agents. The smoothed occupation field $\tilde{n}$ at site $\bm{r}$ is thus
given by the discrete convolution

$\displaystyle\tilde{n}(\bm{r})=\sum_{\bm{r}^{\prime}}G_{\sigma}(\bm{r}-\bm{r}^{\prime})n(\bm{r}^{\prime}).$
(1)

At each time step, an agent $k$ can decide to move out from their occupied
site $\bm{r}_{k}$ and to settle on a new, randomly chosen, empty site
$\bm{r}_{k}^{\prime}$ where the utility $u[\tilde{n}(\bm{r}_{k}^{\prime})]$ –
quantifying the agent’s satisfaction – might exceed their previous utility
$u[\tilde{n}(\bm{r}_{k})]$. We assume that the decision to move to the new
site is a function of the utility difference $\Delta u_{k}\equiv
u[\tilde{n}(\bm{r}_{k}^{\prime})]-u[\tilde{n}(\bm{r}_{k})]$. While the very
existence of the utility function is debatable from a behavioural standpoint
[25], classical economics has traditionally taken agents to be strict utility
maximizers, meaning the move will be accepted if $\Delta u_{k}>0$ and rejected
otherwise. In order to mitigate this assumption, a common approach is to
introduce a stochastic decision rule of the form

$\displaystyle\mathbb{P}(\bm{r}_{k}\to\bm{r}_{k}^{\prime})=f_{\Gamma}(\Delta
u_{k}),$ (2)

where the function $f_{\Gamma}$ is larger than $\frac{1}{2}$ whenever $\Delta
u_{k}>0$. Typically, $f_{\Gamma}$ is a positive and monotonic function of the
utility difference, with $\lim_{x\to-\infty}f_{\Gamma}(x)=0$ and
$\lim_{x\to+\infty}f_{\Gamma}(x)=1$ [26]. The parameter $\Gamma\geq 0$, known
as the intensity of choice, or simply the rationality, quantifies the
propensity of agents to go for strict utility maximizing. In particular,
$\Gamma\to 0$ corresponds to random decision making, while $\Gamma\to\infty$
means perfectly rational agents.

In reality, the specific shape of the function $f_{\Gamma}$ is unknown. In the
socio-economics literature, it is most of the time taken as the logistic
function

$\displaystyle f_{\Gamma}(x)=\frac{1}{1+e^{-\Gamma x}},$ (3)

defining the so-called _logit rule_ [27, 26]. The various reasons and
justifications of this decision rule are discussed and summarized in [21]. In
a nutshell, it can be motivated axiomatically [27], or by the fact that
$f_{\Gamma}$ is a maximum entropy distribution and therefore optimizes an
exploration-exploitation tradeoff when the cost associated with information
scales as $1/\Gamma$ [28, 29]. As empirical evidence supporting this choice
remains extremely scarce, its popularity is in reality largely motivated by
convenience [25]. Indeed, many calculations are made possible thanks to the
fact that it preserves detailed balance with respect to the Gibbs-Boltzmann
measure in the particular case where agents’ utility change also coincides
with a global utility difference [30]. In this context, $T\equiv 1/\Gamma$ can
naturally be interpreted as the temperature, or “social temperature”, of the
system. In the following, the function $f_{\Gamma}$ will be left unspecified,
unless stated otherwise. In the Monte Carlo simulations we will notably use
the logit rule for simplicity.

The last ingredient to specify is the utility function $u$ of the agents. As
stated above, we assume that the utility depends on the locally smoothed
occupation $\tilde{n}$ only, and that it is non monotonic. As in Ref. [15], we
assume that the utility is maximal for some density
$\rho^{\star}\geq\frac{1}{2}$. We specifically choose for the simulations

$\displaystyle u(x)=-\left|x-\rho^{\star}\right|^{\alpha},$ (4)

with $\alpha>0$, see Fig. 1(a), but theoretical computations below will keep
$u$ unspecified.

Figure 1: (a) Utility function $u(\rho)=-|\rho-\rho^{\star}|^{\alpha}$ for
$\rho^{\star}=0.5$, $\alpha=\\{0.5,1,2\\}$. Panels (b), (c) and (d) show
snapshots of the stationary state for these different utility functions,
starting from the same homogeneous profile at $\rho_{0}=0.5$. Here
$\Gamma=100$, $\sigma=3$ and $L_{x}=L_{y}=100$. The stationary density
$\rho_{d}$ in the dense phase is $\rho_{d}=0.575(5)$ for $\alpha=0.5$ in (b),
$\rho_{d}=0.575(5)$ for $\alpha=1$ in (c), and $\rho_{d}=0.585(5)$ for
$\alpha=2$ in (d). These bulk densities are all significantly higher than the
density $\rho^{\star}$ for which agents maximize their utility. Note the
accumulation of agents at the edge of the empty domain in (b) and (c), see
Sec. IV.

###  II.2 In or out of equilibrium?

As mentioned above, an often unspoken motivation for the use of the logit rule
in the modeling of socioeconomic systems is that it may satisfy detailed
balance. Indeed, as described by Grauwin et al. [15, 31] or in [21] in a more
general setting, if one manages to find a system-wide energy-like function
$\mathcal{H}$ such that

$\displaystyle\begin{split}\Delta
u_{k}&=\mathcal{H}([\\{n(\bm{r})\\},n(\bm{r}_{k}^{\prime})=1,n(\bm{r}_{k})=0])\\\
&\quad-\mathcal{H}([\\{n(\bm{r})\\},n(\bm{r}_{k}^{\prime})=0,n(\bm{r}_{k})=1]),\end{split}$
(5)

then the usual tools of equilibrium statistical mechanics can be used. The
steady-state distribution of agents is notably identified as the minimum of
the free energy, which is a Lyapunov function of the dynamics prescribed by
the logit rule.

At the agent level, the existence of such a global quantity is usually the
symptom of either altruistic individuals (that voluntarily maximize some
collective satisfaction) or of a central planner (that constructs individual
rewards towards a collective objective). Outside of these two cases, the
existence of a free energy when agents are individualistic is in fact
restricted to a limited number of carefully chosen models (see [22] for a
related discussion in the context of microeconomics). In the literature of
Schelling-like models, taking a city divided in neighborhoods or blocks [15],
where agents share the same utility, yields such a free energy description
(which is importantly not a simple aggregation of individual utilities). In
our model, however, this is no longer true.

Figure 2: Loop of four configurations with $N=3$ agents on $M=5$ sites
breaking Kolmogorov’s criterion when the utility is non-linear and is a
function of an individual perceived density. Shaded and unshaded nodes
correspond to occupied and empty sites respectively. The dashed line indicates
a possible segmentation of the system into two distinct neighborhoods.

To explicitly show that the dynamics breaks detailed balance despite the logit
rule, one may consider a small system and find a specific cycle breaking
Kolmogorov’s criterion [32]. Such a cycle between four consecutive states with
$N=3$ agents placed on a one-dimensional “city” with $M=5$ sites is
illustrated in Fig. 2. The ratio of transition rates between states $X$ and
$Y$, that differ by an agent located on sites $i$ in $X$, versus $j$ in state
$Y$, is given by

$\frac{W_{X\to Y}}{W_{Y\to
X}}=\frac{1+\mathrm{e}^{-\Gamma[u(\tilde{n}_{i}^{X})-u(\tilde{n}_{j}^{Y})]}}{1+\mathrm{e}^{-\Gamma[u(\tilde{n}_{j}^{Y})-u(\tilde{n}_{i}^{X})]}}=\mathrm{e}^{\Gamma[u(\tilde{n}_{j}^{Y})-u(\tilde{n}_{i}^{X})]}.$
(6)

As a result, the ratio between the product of forward rates, $W_{+}$, and the
product of backwards rates, $W_{-}$, in the cycle shown in Fig. 2, is given by

$\displaystyle\frac{W_{+}}{W_{-}}=\mathrm{e}^{\Gamma[u(\tilde{n}_{5}^{B})-u(\tilde{n}_{3}^{A})-u(\tilde{n}_{2}^{B})+u(\tilde{n}_{3}^{D})+u(\tilde{n}_{2}^{A})-u(\tilde{n}_{5}^{D})]}.$
(7)

For a generic non-linear utility function, $W_{+}\neq W_{-}$, which is a
signature of nonequilibrium dynamics. For a linear utility function on the
other hand, considering that the convolution kernel $G_{\sigma}$ is isotropic,
all terms in the exponential cancel out, leading to $W_{+}=W_{-}$ (which would
be also satisfied for any other cycle). In this situation, the utility
difference can simply be interpreted as an energy difference, where the kernel
$G_{\sigma}$ plays the role of a pairwise interaction potential between the
agents. Interestingly, this small cycle also illustrates how the introduction
of neighborhoods can salvage the equilibrium description for a generic
utility. Splitting the lattice in two neighborhoods along the dashed line
shown in Fig. 2 and taking an identical value of $\tilde{n}$ for all agents on
each neighborhood, the terms in the exponential in Eq. (7) indeed cancel out
for any utility function since $\tilde{n}_{5}^{B}=\tilde{n}_{5}^{D}$,
$\tilde{n}_{3}^{A}=\tilde{n}_{2}^{A}$ and
$\tilde{n}_{2}^{B}=\tilde{n}_{3}^{D}$.

Figure 3: (a) Typical dense domain size $L_{d}(t)$ during coarsening as a
function of time $t$. A unit of time is defined as $N$ Monte Carlo steps,
where $N$ is the number of agents. $L_{d}(t)$ is averaged over 5 independent
simulations. (b), (c), (d) and (e) show snapshots at different times. Starting
from a disordered configuration, we quench the system at low temperature, or
high rationality $\Gamma$, corresponding to $T\simeq T_{c}/6$. Parameters:
$L_{x}=L_{y}=600$, $\rho_{0}=0.3$, $\sigma=1$, $\alpha=3/2$, $T=0.01$.

###  II.3 Microscopic simulations

Having established the out-of-equilibrium nature of our model, we start by
performing numerical simulations to assess whether the concentration of agents
in overly dense regions is generic and robust to different shapes of the
utility function. Here, all numerical simulations are performed on a two-
dimensional grid with periodic boundary conditions. The utility is maximal for
$\rho^{\star}=1/2$. For the sake a simplicity, here we use the logit decision
rule and a truncated Gaussian kernel

$G_{\sigma}(\bm{r})=\begin{cases}\frac{1}{N_{\sigma}}\mathrm{e}^{-\frac{1}{2\sigma^{2}}\|\bm{r}\|^{2}},\quad\text{
if $\|\bm{r}\|$ }\leq 4\sigma,\\\ 0,\quad\text{otherwise},\end{cases}$ (8)

where $N_{\sigma}$ enforces the normalization of the kernel.

####  II.3.1 Phase separation

For large system size $L_{x},L_{y}\gg\sigma$, we explore the behavior for
different global densities $\rho_{0}=N/(L_{x}L_{y})$ and for various
rationality parameters $\Gamma$. Numerical results are qualitatively similar
for all the values of $\alpha$ we tested, ranging from $\alpha=0.5$ to
$\alpha=2$, see Fig. 1. The phenomenology can be summarized as follows. When
rationality is low ($\Gamma\to 0$, $T\to\infty$), the stationary state remains
homogeneous because agents settle at random. When rationality is high, agents
may aggregate in dense clusters, which can surprisingly be more crowded than
what agents’ utilities prescribe. This was already discussed in [15] where the
authors point out that the homogeneous state is actually an unstable Nash
equilibrium, even though all agents maximize their utility. The
destabilization occurs as one agent randomly moves to another region (with no
regard to the effect it may have on the other agents utilities), which
decreases the average density at their original site and increases the average
density where they settle. Agents in the lower-density region will eventually
move to gain the utility they lost when their neighbors moved out. This
dynamics will eventually empty some regions, in which agent’s return becomes
statistically less and less probable. The final state, where a dense phase and
an empty phase coexist, is a stable Nash equilibrium.

One can quantify the condensation dynamics when starting from the homogeneous
state and taking high rationality. The system undergoes a spinodal
decomposition where dense clusters grow and merge until there is one large
dense cluster only, as shown in Fig. 3. The final cluster topology ultimately
depends on noise realization and on the box dimensions. We measure the cluster
size $L_{d}(t)$ as a function of time $t$ using the radial structure factor
(see App. A). We find $L_{d}(t)\sim t^{1/z}$, with the dynamical exponent
$z\in[2,3]$, reminiscent of the coarsening exponent observed in a 2D Ising
system with long-range Kawasaki dynamics [33, 34, 35, 36]. Interestingly, and
consistent with the findings of [36] in the low temperature region, our
results suggest an exponent closer to the local Kawasaki dynamics result $z=3$
(see Fig. 3(a)), despite long-range particle displacements.

####  II.3.2 Critical point and critical exponents

The complete phase separation that occurs when rationality is high indicates
the use of the order parameter $m\equiv\rho_{d}-\rho_{g}$, where
$\rho_{d},\rho_{g}$ are the average densities of the dense and “gas” (dilute)
phases, respectively. At the critical point $(\rho_{c},T_{c})$, we expect a
second-order phase transition where $m$ goes to $0$ with power-law scaling

$\displaystyle m\underset{\tau\to 0^{+}}{\sim}\tau^{\beta},$ (9)

where $\tau=(T_{c}-T)/T_{c}>0$ defines the rescaled temperature difference,
and $\beta$ is the order-parameter critical exponent. Measuring the critical
exponents allows one to determine to which universality class the system
belongs to, providing precious information on the system behavior at large
scales. Since simulations are carried out in finite systems, measuring the
critical point with precision requires numerical ruse. We follow the approach
that has been extensively used to measure critical exponents in systems
undergoing a Motility-Induced Phase Separation (referred to as MIPS) [37, 38,
39], see App. B.

Simulations are performed in a rectangular domain of size $L_{x}\times L_{y}$,
with $L_{x}=3L_{y}$, with periodic boundary conditions to keep flat interfaces
between a stripe of liquid (dense phase) and a stripe of gas (dilute phase).
Starting with the dense phase in the center of the system, we track the center
of mass such that we always compute the densities in the bulk of each phase.
To compute the local density inside the bulk of each phase, we consider square
boxes of size $\ell=L_{y}/2$, centered either in $0$ in the gas bulk or
centered in $L_{x}/2$ in the dense bulk (Fig. 8). The local density in each
box fluctuates and it is given by $\rho_{b}=N_{b}/\ell^{2}$ with $N_{b}$ the
number of agents in the box $b$ in a given realization of the system. The
distribution of the density in the system is thus bimodal for $T<T_{c}$ and
unimodal when the system is homogeneous. Defining

$\displaystyle\Delta\rho=\frac{N_{b}-\langle N_{b}\rangle}{\ell^{2}},$ (10)

where the $\langle\cdot\rangle$ stands for averaging on the four boxes and on
independent realizations of the simulation, we compute the celebrated Binder
cumulant [40, *rovere_simulation_1993, 42]

$\displaystyle
Q_{\ell}(\Delta\rho,T)=\frac{\langle(\Delta\rho)^{2}\rangle}{\langle(\Delta\rho)^{4}\rangle},$
(11)

for a given box size $\ell$ and a given temperature $T$. For $\ell$ large
enough, the curves $Q_{\ell}(T)$ all intersect in $T=T_{c}$ where the behavior
of the system is universal. It is important to mention that the critical
density is not known _a priori_. It has to be assessed beforehand to ensure
that the system, as $T$ changes, goes through the critical point, where the
phase transition is of second-order type. To locate $\rho_{c}$, we compute the
Binder cumulant at fixed temperature, close to the estimated critical point,
for various densities $\rho_{0}$. The critical density then corresponds to the
maximal fluctuations of $\Delta\rho$, translated in a peak of the Binder
cumulant, see Fig. 4(a). Once the critical point is precisely located,
additional critical exponents can be measured. Notably, defining the
susceptibility $\chi$ as

$\displaystyle\chi\equiv\frac{\langle(N_{b}-\langle
N_{b}\rangle)^{2}\rangle}{\langle
N_{b}\rangle}=\frac{\langle(\Delta\rho)^{2}\rangle}{\langle
N_{b}\rangle}\ell^{4},$ (12)

one obtains

$\displaystyle\chi\sim\ell^{\gamma/\nu},\quad\frac{dQ_{\ell}}{d\tau}\Big{|}_{\tau=0}\sim\ell^{1/\nu},$
(13)

at the critical point.

Figure 4: Numerical experiments for $\sigma=1$, $\alpha=3/2$. (a) Binodal
densities measured for $L_{x}=200$ and $L_{y}=66$ ($\ell=33$), inset showing
the Binder cumulant as a function of the density and fitted (continuous line)
to determine the critical density. (b), (c) and (d) show the numerical
measurements of the critical exponents close to the critical point
$(\rho_{c},T_{c})=(0.271,0.0620)$ determined using various system sizes
ranging from $\ell=20$ to $\ell=40$.

We report in Fig. 4 the various results on the critical point and on the
critical exponent for $\sigma=1$ and $\alpha=3/2$. Using the Binder cumulant,
one identifies the critical point at $\rho_{c}=0.271(5)$ and
$T_{c}=0.0620(2)$, where the uncertainty on the last digit appears in the
parentheses. The phase diagram in space $(\rho,T)$ is shown in Fig. 4(a), the
black star indicates the critical point and the circular markers show the
densities of the coexisting phases: they define the _binodal_ frontier. The
exponent $\beta$ is directly measured from the order parameter $m$ as function
of reduced temperature $\tau$, at a fixed system size $L_{x}=220$ [see
Fig.9(c)]. From the Ising-2D ansatz, we check that $\nu=1$ yields a neat
collapse of the Binder cumulant, see Fig. 9(b). The exponent $\gamma$ is
obtained by varying the system size at the critical temperature $T_{c}$ and
assuming $\nu=1$ [see Fig.9(d)]. We report in Table 1 the values found for the
critical exponents in the cases $\alpha=3/2$ (Fig. 4) and $\alpha=1/2$ (not
shown here).

Model | $\rho_{c}$ | $T_{c}$ | $\beta$ | $\gamma$
---|---|---|---|---
Ising 2D (exact) | 0.5 |  | 0.125 | 1.75
$\alpha=1/2$, $\sigma=1$ | 0.309(5) | 0.0983(5) | 0.120(8) | 1.71(5)
$\alpha=3/2$, $\sigma=1$ | 0.271(5) | 0.0620(2) | 0.119(5) | 1.74(5)
Table 1: Critical density and exponents for nonequilibrium Sakoda-Schelling
model for $\alpha=1/2$ and $\alpha=3/2$.

They differ by less than 5% from the 2D Ising static exponents. These results
enjoin us to assert with a high degree of confidence that the model considered
here belongs to the 2D-Ising universality class. Since the system is out of
equilibrium and particle displacements can be of infinite range, recovering
the Ising universality class is a priori nontrivial. However, finding other
critical exponents would have been surprising since the ingredients at play
are the ones of the Ising model, namely, short-range and isotropic
interactions, a homogeneous medium and a two state degree of freedom (sites
are empty or occupied). This result must also be put into perspective with the
recent debate on the universality class(es) of systems undergoing MIPS [37,
39, 38, 43, 44], and their associated active field theories [45, 46]. Notably
here, our interaction kernel $G_{\sigma}$ provides a so-called _quorum-
sensing_ interaction, like that found in assemblies of bacteria [47]. The
particle dynamics is however quite different for bacteria and for our agents.
The remaining of the paper shall be devoted to establishing a quantitative
relation between our Sakoda-Schelling occupation model and the field-theory
descriptions of MIPS. The first step along this path is to formulate a mean-
field approximation of our model.

##  III Field theory and the local-move approximation

###  III.1 General description

The computation starts by writing the expectation of the occupation number
$n_{\bm{r},s+1}\equiv n(\bm{r},s+1)$ of site $\bm{r}$ at time $s+1$,
conditioned on the previous configuration $\\{n_{\bm{r},s}\\}$. Averaging over
multiple realizations of noise and using a mean-field approximation in which
all correlation functions factorize, one obtains

$\displaystyle\begin{split}\langle n_{\bm{r},s+1}\rangle-\langle
n_{\bm{r},s}\rangle&=(1-\langle
n_{\bm{r},s}\rangle)\sum_{\bm{r}^{\prime}\neq\bm{r}}\langle
n_{\bm{r}^{\prime},s}\rangle f_{\Gamma}(\Delta
u^{s}_{\bm{r}^{\prime}\to\bm{r}})\\\ &-\langle
n_{\bm{r},s}\rangle\sum_{\bm{r}^{\prime}\neq\bm{r}}(1-\langle
n_{\bm{r}^{\prime},s}\rangle)f_{\Gamma}(\Delta
u^{s}_{\bm{r}\to\bm{r}^{\prime}}),\end{split}$ (14)

where $\Delta u^{s}_{\bm{r}\to\bm{r}^{\prime}}\equiv
u(\langle\tilde{n}_{\bm{r}^{\prime},s}\rangle)-u(\langle\tilde{n}_{\bm{r},s}\rangle)$.
For convenience, we take the continuous time and continuous space limit,
following the common procedure to obtain a mean-field description of exclusion
processes on lattices (see e.g. [48]). The average occupation number $\langle
n\rangle$ is now described by the density $\rho$, while the spatially smoothed
average occupation number $\langle\tilde{n}\rangle$ is described by the field
$\phi\equiv G_{\sigma}*\rho$. The master equation for the occupation
probability then takes the form of a noiseless hydrodynamic equation, in our
case:

$\displaystyle\partial_{t}\rho(x,t)=$
$\displaystyle[1-\rho(x,t)]\int\mathrm{d}y\,\rho(y,t)w_{\Gamma}([\phi],y,x,t)$
(15)
$\displaystyle-\rho(x,t)\int\mathrm{d}y\,[1-\rho(y,t)]w_{\Gamma}([\phi],x,y,t),$

with the transition rate from $y$ to $x$ explicitly given by

$\displaystyle w_{\Gamma}([\phi],y,x,t)=\omega
f_{\Gamma}\left(u(\phi(x,t))-u(\phi(y,t))\right),$ (16)

where $\omega$ is homogeneous to an inverse time scale and $f_{\Gamma}$ is
left unspecified. Equation (15) is valid in any dimension, but, for
simplicity, we will work out the mean-field computations in dimension 1 in
space. This can be justified _a posteriori_ when we compare the mean-field
(MF) to the Monte Carlo (MC) simulations. Let us also mention that the
dimension does not play a role in determining the phase densities in the
steady state of coarse-grained field theories (Allen-Cahn [49], Cahn-Hilliard
[50], etc.).

Integrating Eq. (15) over space, one immediately sees that the total density
$\int\rho$ is conserved. One can also check that in the very specific case
where $u(\phi)$ is linear in $\phi$, one can build a free-energy functional
that is a Lyapunov function of the non-local MF dynamics, ensuring a
convergence towards local minima and preventing limit cycles and oscillatory
dynamics. This is a natural consequence of the fact that detailed balance is
satisfied at microscopic level. In App. C, we construct this free energy and
show that the dynamics is relaxational.

###  III.2 Linear stability analysis

In the general case, we would like to understand how the homogeneous state
becomes unstable. To do so, we consider a small perturbation around the
homogeneous state: $\rho(x,t)=\rho_{0}+\rho_{1}(x,t)$, with $\rho_{1}$ the
perturbation. By linearity of the convolution, one has
$\phi(x,t)=\rho_{0}+\phi_{1}(x,t)$, with $\phi_{1}\equiv G_{\sigma}*\rho_{1}$.
A Taylor expansion of Eq. (15) combined with mass conservation (i.e
$\int_{D}\rho_{1}=\int_{D}\phi_{1}=0$, where $D$ is the full domain), finally
yields:

$\displaystyle\begin{split}\partial_{t}\rho_{1}(x,t)=&\
2\Omega\rho_{0}(1-\rho_{0})f^{\prime}_{\Gamma}(0)u^{\prime}(\rho_{0})\phi_{1}(x,t)\\\
&-\Omega f_{\Gamma}(0)\rho_{1}(x,t),\end{split}$ (17)

with $\Omega$ the full domain size. Defining the Fourier transform for any
field $h$ as $\hat{h}(k)=\int\mathrm{d}x\,e^{-ikx}h(x)$, one obtains

$\displaystyle\partial_{t}\hat{\rho}_{1}(k,t)=\Lambda(k)\hat{\rho}_{1}(k,t),$
(18) $\displaystyle\Lambda(k)=\Omega
f_{\Gamma}(0)\left(2\rho_{0}(1-\rho_{0})\frac{f^{\prime}_{\Gamma}(0)}{f_{\Gamma}(0)}u^{\prime}(\rho_{0})\hat{G}_{\sigma}(k)-1\right).$
(19)

This last equation shows that the homogeneous state is unstable if there
exists a mode $k^{\star}$ such that

$\displaystyle
2\rho_{0}(1-\rho_{0})\frac{f^{\prime}_{\Gamma}(0)}{f_{\Gamma}(0)}u^{\prime}(\rho_{0})\hat{G}_{\sigma}(k^{\star})>1.$
(20)

The manifold for which the inequality becomes an equality defines the spinodal
in the phase diagram $(\rho_{0},\Gamma)$. In particular, for any monotonically
decreasing kernel $G_{\sigma}(|x|)\in L_{2}(\mathbb{R})$, one has
$\hat{G}_{\sigma}(0)>|\hat{G}_{\sigma}(k)|$, such that for large system size,
the stability of the homogeneous state is given by the stability of modes
$k\to 0$, and the spinodal is thus defined by the equation

$\displaystyle
2\rho_{0}(1-\rho_{0})\frac{f^{\prime}_{\Gamma}(0)}{f_{\Gamma}(0)}u^{\prime}(\rho_{0})=1.$
(21)

Note that this criterion is generic as it only depends on the decision rule
through $f_{\Gamma}(0)$ and $f^{\prime}_{\Gamma}(0)$. The simulations also
reveal the existence of a bistable region in the vicinity of this spinodal.
This is the binodal region, where hysteresis and bistability can notably be
observed, and which can be fully characterized in the case of an equilibrium
system [12]. Here however, there is a priori no free energy one can rely on to
describe the nucleation scenario and to obtain the densities of the phase-
separated state.

Figure 5: Comparison between Monte Carlo simulations and mean-field results
for $\alpha=3/2$, $\sigma=7$ and $L_{x}=200$, $L_{y}=66$ ($\ell=33$). (a)
Phase diagram in the $(\rho_{0},T)$ plane. The mean-field binodal (continuous
black line) is given by measuring the densities of the bulk of each plateau in
a phase separated state. The circles are the bulk averaged densities in Monte
Carlo (MC) simulations. The dashed black line represent the mean-field
spinodal, which is obtained analytically from the linear stability analysis
(see Eq. (20)), with $\hat{G}_{\sigma}(k_{1})=e^{-\sigma^{2}k_{1}^{2}/2}$ and
$k_{1}=2\pi/L_{x}$. The diamonds indicate the loss of stability of the
homogeneous state in the MC simulations. The black square is the critical
point for $\sigma/L\to 0$. (b) Averaged density profile $\rho(x)$ from MC
simulations (continuous green line) for $\rho_{0}=0.35$, $T=0.05$. The dashed
black line is the stationary solution of the mean-field equation Eq. (15) for
the same parameters, solved on a grid of step size 1 with a Euler explicit
scheme.

###  III.3 Comparison to microscopic simulations

The MF prediction is expected to be accurate for systems with high
connectivity, which here corresponds to large $\sigma$. In the following, we
shall take the limit $L\to\infty$, $\sigma\to\infty$ with $\sigma/L\to 0$ to
obtain mean-field predictions that are independent of both $\sigma$ and $L$,
and perform numerical simulations as close as possible to this scaling regime.

The first analytical prediction of the MF description is the spinodal, that
determines the onset of instability of the homogeneous state, see Eq. (21).
The spinodal is the dashed line in the $(\rho,T)$ phase diagram in Fig. 5(a).
To check the prediction, we start in the MC simulations from a uniformly
distributed configuration of agents for three different values of temperature,
$T=0.04$, $0.08$, $0.11$, and we detect the frontier across which the
homogeneous profile either coarsens, or needs a nucleation event to converge
to the separated state. This frontier is marked with the diamonds, which
agrees with the MF prediction.

Second, the MF dynamical Eq. (15) can be solved numerically with an Euler
explicit scheme. From the numerical solution, one obtains the densities of the
bulk of each phase when a phase separation occurs: these densities define the
binodal, the continuous line in Fig. 5(a). These MF phase densities are
perfectly recovered by the MC simulations (circles). In addition, one can
compare the steady-state average density profile from MC simulations to the
mean-field stationary density, which superimpose almost exactly, see Fig.
5(b).

As previously stated, the MF predictions fail for small values of $\sigma$.
The phase diagram in Fig. 4(a) is for instance obtained for $\sigma=1$, and
indeed strongly differs from the MF solution. For $\sigma=1$, we notably
identify the critical point at $(\rho_{c},T_{c})=(0.271,0.0620)$, whereas the
MF predicts $(\rho_{c},T_{c})_{\mathrm{MF}}=(0.2763,0.1418)$, where, as
expected, $T_{c}^{\sigma=1}<T_{c}^{\mathrm{MF}}$.

###  III.4 Local-move approximation

To make progress into the identification of a possible effective free energy
functional, it may be convenient to consider slightly modified dynamics where
jumps are now only authorized in the direct neighborhood of the agents.
Indeed, considering an evolution enforcing a local mass conservation will
allow for more familiar partial differential equations (PDEs) and field
theoretic approaches on conserved scalar fields. Here, the absence of
macroscopic density currents in the steady state, both in MC simulations and
in the MF solution suggests that the system generically converges to a
stationary stable fixed point, where the details of the dynamics become
inconsequential. In addition, when the majority of agents have aggregated in a
single dense cluster in the steady state, it is unlikely that they would
perform moves outside of the bulk, in low-density regions, since the utility
there is minimal. The local-move approximation, as it strongly simplifies the
description, thus appears natural.222Dynamically, the coarsening exponent
$z\simeq 3$ displayed in Fig. 3(a), and which is also observed in a Cahn-
Hilliard relaxation dynamics can also be invoked to support the idea of local
moves.

Following the Taylor expansion outlined in App. D, the local mean-field
dynamics is given by

$\displaystyle\partial_{t}\rho$
$\displaystyle=f_{\Gamma}(0)\partial_{x}^{2}\rho-2f_{\Gamma}^{\prime}(0)\partial_{x}[\rho(1-\rho)\partial_{x}u],$
(22)

which can be rewritten as the canonical equation for the mass-conserving
dynamics

$\partial_{t}\rho=\partial_{x}[M[\rho]\partial_{x}\mu([\rho],x)],$ (23)

with the mobility operator $M[\rho]=\rho(1-\rho)$, stemming from the non-
overlaping nature of the agents, and with the chemical potential
$\mu=\mu_{\mathrm{ent.}}+\mu_{\mathrm{util.}}$ where

$\displaystyle\mu_{\mathrm{ent.}}$
$\displaystyle=f_{\Gamma}(0)\log\left(\frac{\rho}{1-\rho}\right)$ (24)
$\displaystyle\mu_{\mathrm{util.}}$
$\displaystyle=-2f_{\Gamma}^{\prime}(0)u[\phi(x)].$ (25)

The first contribution to the chemical potential $\mu_{\mathrm{ent.}}$ is
purely local and accounts for entropy in the system where agents cannot
overlap. The second contribution $\mu_{\mathrm{util.}}$ encodes the drive from
agents’ utility. This term exhibits non-locality with respect to the field
$\rho$, and as a consequence, cannot be expressed as a functional derivative
of any free energy, in general [51, 52, 53]. However, in the particular case
of a linear utility in $\phi$, one again recovers that
$\mu_{\mathrm{util.}}+\mu_{\mathrm{ent.}}$ can be written as the functional
derivative of the free energy $\mathcal{F}$ given in App. C and, as a
consequence, the dynamics (23) becomes a gradient descent [54]. Let us
emphasize that, here again, the decision rule is kept general, and that the
entire local dynamics only depend on it through $f_{\Gamma}(0)$ and
$f^{\prime}_{\Gamma}(0)$.

Performing the linear stability analysis on the dynamics with local moves (see
App. D), we find that the criterion for the homogeneous solution to be
unstable is identical to that given in Eq. (21), when moves are global. Also,
the stationary density profiles computed either with the local, or with the
non-local MF PDEs for the same parameters are identical, as shown in App. E.
Both these observations therefore allow us to confirm the relevance of the
local-move approximation to characterize the system in the long-time limit.

Finally, note that the local hydrodynamic equations can also be obtained using
the path integral approach on a lattice [55], which, in passing, provides the
fluctuating hydrodynamics:

$\displaystyle\partial_{t}\rho=\partial_{x}\left[\rho(1-\rho)\partial_{x}\mu([\rho],x)+\sqrt{\rho(1-\rho)}\xi\right],$
(26)

where $\xi(x,t)$ is a Gaussian white noise with zero mean and with
$\langle\xi(x,t)\xi(x^{\prime},t^{\prime})\rangle=2f_{\Gamma}(0)\delta(t-t^{\prime})\delta(x-x^{\prime})$.
We then remark that when the utility is linear, the stochastic field evolution
describes a complete equilibrium dynamics, irrespective of the choice of the
decision rule: A rule that breaks detailed balance at the microscopic level
can still lead to an equilibrium field theory after coarse graining. Similar
findings had been pinpointed in active matter models [56, 52]. While not
central to the present work, the fluctuations can be studied in more detail,
providing information on the nucleation scenarii and on transition paths
between macroscopic states for instance [57, 58, 59, 60]. The study of the
associated functional Fokker-Planck equation using the tools described in [61,
62] may also be an interesting perspective for future works. In the case of
non-local moves, the formalism from [55] cannot be straightforwardly adapted,
since the local gradient expansion of the jump rates in the action breaks
down. Establishing an appropriate fluctuating hydrodynamic description in the
case of non-local dynamics is therefore an open problem.

##  IV Generalized-thermodynamic construction

Since the previous section has shown that the phase separation is well
described by the local-move approximation, we can now use the machinery of
field theory for scalar active matter (e.g. Active Model B), as developed in
[23, 63, 64]. This mapping will notably allow us to obtain the binodal
densities for some class of utility functions detailed below.

###  IV.1 The generalized-thermodynamic expansion

Even though $\mu$ in Eq. (23) cannot be written as the functional derivative
[52, 51, 53], the dynamics can be analyzed by resorting to a gradient
expansion. Indeed, expanding the chemical potential up to
$O(\nabla^{4},\rho^{2})$ terms yields

$\displaystyle\mu[\rho]=g_{0}(\rho)+\lambda(\rho)(\nabla\rho)^{2}-\kappa(\rho)\nabla^{2}\rho+O(\nabla^{4},\rho^{2}),$
(27)

with $g_{0}$, $\lambda$, $\kappa$ local function of the field $\rho$, and a
generalized thermodynamic mapping [63, 65] can yield the prediction of the
binodal densities.

For simplicity, we will now assume that the smoothing kernel is a Gaussian
distribution of zero mean and variance $\sigma^{2}$. In Fourier space, the
smoothed field is given by
$\hat{\phi}(k)=\hat{\rho}_{k}\exp({-{\sigma^{2}k^{2}}/{2}})$, which can be
truncated to leading order:

$\displaystyle\hat{\phi}_{k}\simeq\hat{\rho}_{k}\left(1-\frac{k^{2}\sigma^{2}}{2}+O(\sigma^{4}|k|^{4})\right).$
(28)

In real space, this translates into
$\phi=\rho+\frac{\sigma^{2}}{2}\nabla^{2}\rho+O(\nabla^{4},\rho)$. This allows
us to further expand the $\mu_{\mathrm{util.}}$ given in Eq. (25). To leading
order in the $O(\nabla,\rho)$ expansion, one has

$\displaystyle\mu_{\mathrm{util.}}=-2f^{\prime}_{\Gamma}(0)\left[u(\rho)+\frac{\sigma^{2}}{2}u^{\prime}(\rho)\partial_{x}^{2}\rho+O(\partial_{x}^{4},\rho)\right].$
(29)

Combining this expansion of $\mu_{\mathrm{util.}}$ with the entropic
contribution $\mu_{\mathrm{ent.}}$, it is now possible to identify the
different terms in Eq. (27), namely:

$\displaystyle
g_{0}(\rho)=-2f^{\prime}_{\Gamma}(0)u(\rho)+f_{\Gamma}(0)\log\left(\frac{\rho}{1-\rho}\right);$
(30)
$\displaystyle\lambda(\rho)=0;\quad\kappa(\rho)=f^{\prime}_{\Gamma}(0)\sigma^{2}u^{\prime}(\rho).$
(31)

This identification enables us to follow up to the next step, which is finding
the proper function $R(\rho)$ and the generalized functional $\mathcal{G}[R]$
by means of which the dynamics will be given by

$\displaystyle\partial_{t}\rho(x,t)=\partial_{x}\cdot\left[M[\rho]\partial_{x}\frac{\delta\mathcal{G}}{\delta
R(x,t)}\Big{|}_{R(\rho)}\right].$ (32)

A double-tangent construction on $\mathcal{G}[R]$ then provides the binodal
densities [63]. Since $\lambda(\rho)=0$, the differential equation that the
function $R$ must satisfy (see [63, 65]) is

$\displaystyle\kappa(\rho)R^{\prime\prime}(\rho)=-\kappa^{\prime}(\rho)R^{\prime}(\rho),$
(33)

which simplifies into $(\kappa R^{\prime})^{\prime}=0$, where the ′ denotes
the derivative with respect to $\rho$.

###  IV.2 Implications of non-monotonous utilities

The previous equation suggests $R^{\prime}(\rho)=C/\kappa(\rho)$, with $C$
some constant. However, one has to be careful at this stage. In the case
considered here, where the utility of agents reaches its maximum for some
density $\rho^{\star}$, it is clear that $\kappa(\rho)$ undergoes a sign
change at $\rho^{\star}$, and more precisely since $f^{\prime}_{\Gamma}(0)>0$,
we have $\mathrm{sign}[\kappa(\rho)]=\mathrm{sign}[u^{\prime}(\rho)]$. To our
knowledge, the fact that $\kappa(\rho)$ may not remain strictly positive has
never been considered in the active matter literature, even though it bears
important physical meaning. Consider a system of quorum-sensing moving
bacteria whose microscopic velocity $v(\rho)$ is density dependent [47].
Coarse-graining typically yields $\kappa(\rho)\simeq-v^{\prime}(\rho)/v(\rho)$
[66, 65], and one obtains $\kappa(\rho)>0$ when the velocity of the particles
is a decreasing function of local density, eventually leading to bacteria
condensation, i.e. MIPS [56, 67]. A positive $\kappa>0$ is thus naturally
interpreted as an “effective surface tension” in this framework. On the other
hand, a negative $\kappa(\rho)$ would be the reflect of an increasing motility
with increasing bacterial density, which is also biologically relevant if one
considers that bacteria are likely to avoid competition for resources in
crowded areas. Yet, and this is a key remark here, it does not necessarily
mean that the phase separation is arrested or that the system undergoes a
microphase separation when $\kappa<0$, notably because higher order gradient
terms that were discarded in the field expansion then become relevant and may
stabilize the interfaces. More specifically here, for $\rho>\rho^{\star}$,
$u^{\prime}(\rho)<0$ such that the term $O(\partial_{x}^{2}\rho)$ destabilizes
the interfaces but the term $O(\partial_{x}^{4}\rho)$ prevents the gradient
from diverging, with a role similar to the positive bending stiffness in the
Helfrich Hamiltonian of membranes [68, 69]. This leads to density overshoots
and spatial oscillations at the edges of the densely populated domain, but the
dense phase ultimately reaches a plateau, as can be observed in snapshots
(Figs. 1, 5 and 10). Conversely, a strictly positive $\kappa(\rho)$ suppresses
the density overshoots, as illustrated in Fig. 11(b). More generally, we
believe that such a macroscopic feature, namely, a spatially oscillating
density profile could interestingly be exploited in experimental systems to
provide important clues on the microscopic properties of the constituents
under study.

Coming back to $(\kappa R^{\prime})^{\prime}=0$ with $\kappa(\rho^{\star})=0$
and assuming that $R(\rho)$ simply needs to be bijective and continuous, one
deduces that different constants are now a priori needed on each interval
$(0,{\rho^{\star}})$ and $({\rho^{\star}},1)$ to compute $R(\rho)$. For the
specific class of utility function $u(\rho)=-|\rho-{\rho^{\star}}|^{\alpha}$,
$\alpha\neq 2$, one obtains

$\displaystyle
R^{\prime}(\rho)=\begin{cases}\frac{C_{1}}{\sigma^{2}f^{\prime}_{\Gamma}(0)\alpha({\rho^{\star}}-\rho)^{(\alpha-1)}}\text{
if $\rho<{\rho^{\star}}$,}\\\
-\frac{C_{2}}{\sigma^{2}f^{\prime}_{\Gamma}(0)\alpha(\rho-{\rho^{\star}})^{(\alpha-1)}}\text{
if $\rho>{\rho^{\star}}$}.\end{cases}$ (34)

For $\alpha\neq 2$, the function $R$ is given by

$\displaystyle
R(\rho)=\begin{cases}\displaystyle\frac{C_{1}({\rho^{\star}}-\rho)^{2-\alpha}}{\sigma^{2}f^{\prime}_{\Gamma}(0)\alpha(\alpha-2)}+C_{3}\text{
if $\rho<{\rho^{\star}}$,}\\\
\displaystyle\frac{C_{2}(\rho-{\rho^{\star}})^{2-\alpha}}{\sigma^{2}f^{\prime}_{\Gamma}(0)\alpha(\alpha-2)}+C_{4}\text{
if $\rho>{\rho^{\star}}$},\end{cases}$ (35)

with $C_{i}$ the interval dependent constants. These equations show that the
case $\alpha<2$ may admit an acceptable change of variable. For $\alpha\geq
2$, the function $R$ displays a divergence at $\rho={\rho^{\star}}$, which
makes impossible to recover an homeomorphism $R:\rho\mapsto R(\rho)$ on the
whole domain $(0,1)$.

In the next paragraph, we detail the procedure introduced in [63, 65] to
recover the binodal densities. We show that we can extend the procedure to
negative $\kappa(\rho)$ – assuming that higher order gradient terms stabilize
the interface. This is one of the central results of this paper.

The function $R$ is bijective, and can thus be inverted to get $\rho(R)$, and
this allows us to define a new chemical potential $g[R]\equiv\mu[\rho(R)]$.
The functional $\mathcal{G}[R]$ is then obtained by integrating $g[R]$ on each
domain and by gluing together the two integrated parts at ${\rho^{\star}}$.
Explicitly using the notations introduced in [65], we have

$\displaystyle g=\frac{\delta\mathcal{G}}{\delta
R},\qquad\mathcal{G}=\int\mathrm{d}x\left[\Phi(R)+\frac{\kappa}{2R^{\prime}}(\partial_{x}R)^{2}\right],$
(36)

where $\Phi(R)$ defines a generalized free energy density verifying

$\displaystyle\frac{\mathrm{d}\Phi}{\mathrm{d}R}=g_{0}[\rho(R)],$ (37)

with $g_{0}$ defined in Eq. (30). The double-tangent construction on $\Phi(R)$
then yields the binodal densities. To be more explicit, from Eq. (35), we
assume that $R$ is simply given by

$\displaystyle R(\rho)=\begin{cases}-(\rho^{\star}-\rho)^{2-\alpha}\text{ if
$\rho<{\rho^{\star}}$,}\\\ (\rho-\rho^{\star})^{2-\alpha}\text{ if
$\rho>{\rho^{\star}}$},\end{cases}$ (38)

which can be inverted into
$\rho(R)=\rho^{\star}+\operatorname*{sgn}(R)|R|^{\frac{1}{2-\alpha}}$.
Injecting $\rho(R)$ in Eq. (30), one obtains

$\displaystyle\begin{split}g_{0}(R)=&\,2f^{\prime}_{\Gamma}(0)|R|^{\frac{\alpha}{2-\alpha}}\\\
&+f_{\Gamma}(0)\log\left[\frac{\rho^{\star}+\operatorname*{sgn}(R)|R|^{\frac{1}{2-\alpha}}}{1-\rho^{\star}-\operatorname*{sgn}(R)|R|^{\frac{1}{2-\alpha}}}\right].\end{split}$
(39)

The explicit formula for $\Phi(R)$ is more involved and we choose to display
$R(\rho)$ in Fig. 6(a) for $\alpha=3/2$ and $\rho^{\star}=1/2$. To obtain the
binodal densities, one either performs the double-tangent construction on
$\Phi$ (Fig. 6(b), inset) or the Maxwell construction, the latter being easier
to handle numerically [65]. We provide details on these two constructions in
Appendix G. In a nutshell, one obtains the coexistence densities of the dilute
and the dense phases from the constraints of the steady state. Indeed, in the
steady state, the interface between the phases does not move and each phase
has a fixed density, which translates into equality across phases of pressure
and chemical potential, respectively. We then compare the predictions to the
phase densities measured in MC simulations in Fig. 6(b), and report an
excellent match. Note that the interaction range $\sigma$, that appears in
$\kappa$ but not in $\Phi$, does not play a role in the predicted coexistence
densities of the infinite size system, confirming that the sub-optimal
aggregation of agents in a dense cluster is not limited to finite size
lattices.

Figure 6: (a) Utility $u(\rho)$ (solid line), and change of variable $R(\rho)$
(dashed line) for $\alpha=3/2$ and $\rho^{\star}=1/2$. (b) Comparison between
the semi-analytical prediction (dark line) and the binodal densities both
obtained via the Monte Carlo simulations (green circles) and solving
numerically the mean-field Eq. (15) (green dot-dashed line). The decision rule
here is the logit function, whose values in $0$ are $f_{\Gamma}(0)=1/2$ and
$f^{\prime}_{\Gamma}(0)=\Gamma/4$. Inset: double-tangent construction on
$\Phi(R)$ for $T=0.083$. The dense phase is given by $R_{d}=0.170$, yielding
$\rho_{d}=0.53$, and the gaseous phase is given by $R_{g}=-0.69$, yielding
$\rho_{g}=0.022$.

As a final remark, and to provide some validity criterion for the extended
change variable, we should mention that $R(\rho)$ taken alone does not contain
information on the sign of $\kappa(\rho)$. This is why we claim that, if the
system indeed undergoes a true phase separation and that $\nabla\rho$ has
finite left and right limits at each point of the domain, then the sign of
$\kappa$ does not matter and a computation with negative $\kappa$ still
predicts the correct binodal densities. In other words, it is not the negative
sign of $\kappa$ _per se_ that leads to the failure of the gradient expansion,
but rather the fact that there exists some points $x_{c}$ in the domain such
that $|\partial_{x}\rho\small|_{x_{c}^{\pm}}|=+\infty$. The fact that
$\partial_{x}\rho$ may not be continuous but keeps finite values at the right
and at the left of each point of the domain does not seem to be an issue when
performing the gradient expansion, probably because the set of points where
the derivative is higher than some arbitrary value is a set of zero measure.
The breakdown of the gradient expansion for a positive and monotonic
$\kappa(\rho)$ is illustrated in App. F.

##  V Further socioeconomic considerations

###  V.1 Two populations

A natural extension of the problem is to restore some diversity among agents,
as initially considered by both Sakoda and Schelling. Here we consider two
types of interacting agents (say $A$ and $B$), with possibly different utility
functions, which could for example represent higher and lower revenue
individuals, or city dwellers and business professionals, etc.333The recent
work [70] on urban segregation in the United States has brought to our
attention the existence of surveys [71, 72, 73] confirming the idea that
different sub-populations may require markedly different utility functions. A
central question in this case is whether the system reaches fixed points, or
if more complicated dynamics can persist in the long time limit, especially if
the two populations have competing interests. Recent work has been devoted to
studying nonreciprocal interactions between different kinds of particles,
exhibiting the wealth of possible dynamical behavior when particle
displacements are local [74, 75]. An interesting question in our setup is for
instance: do propagating waves (or frustrated states) survive when nonlocal
moves are allowed? Indeed, one may expect that enforcing local displacement
constitutes a dynamical constraint that drives the system in a particular way.
Allowing for nonlocal moves may change the dynamics of how the frustrated
states are resolved. One may think of three major types of interactions:

  * •

First, a cooperative interaction where agents $A$ and agents $B$ may maximize
their utility when agents of opposite type are found in their neighborhood.
This kind of interaction will typically lead to homogeneous well-mixed
systems, or to some condensation into a dense phase where agents are well-
mixed, but since frustration is not implemented in the microscopic rules, we
reasonably expect stationary states.

  * •

Second, each agent type may decide to settle among peers and/or avoid agents
of the other type in their surroundings. One should then expect a complete
phase separation into two domains, one displaying a majority of $A$s and, the
other, a majority of $B$s. Whether the $A-B$ phase separation additionally
displays some condensation depends on the self-affinity of each agent type.

  * •

Third, frustrated situations in which $A$ settles with $A$ but wants to avoid
$B$ agents, while $B$ agents would like to gather and settle close to $A$. In
this situation, we may expect non stationary patterns, stemming from the fact
that all agents cannot be satisfied at the same time.

With this last situation in mind, we have considered the following utility
functions ($u_{A}$ for $A$ agents and $u_{B}$ for $B$ agents):

$\displaystyle u_{A}(x,[\phi_{A,B}])$
$\displaystyle=-|\phi_{A}(x)-\rho^{\star}|^{2}+c_{1}\phi_{B}(x)$ (40)
$\displaystyle u_{B}(x,[\phi_{A,B}])$
$\displaystyle=-|\phi_{A}(x)-\rho^{\star}|^{2}+c_{2}\phi_{B}(x),$ (41)

where $c_{1}<0$ translates the fact that $A$s are fleeing from $B$, and
$c_{2}>0$ translates the fact that $B$s have a tendency to gather with $B$s.
The term $-|\phi_{A}-\rho^{\star}|^{2}$ enjoins both populations to settle
among $A$ populated areas. Of course, the specific shape of utilities taken
here may be restrictive and can be easily generalized.

The extension of the mean-field dynamics to this two population problem is
rather straightforward. Writing $\rho_{A}(x,t)$ (resp. $\rho_{B}(x,t)$) the
density of agents $A$ (resp. $B$) at location $x$ and time $t$, and denoting
the total density by $\rho(x,t)\equiv\rho_{A}(x,t)+\rho_{B}(x,t)$, we now have
an evolution equation of the form

$\displaystyle\partial_{t}\rho_{A}(x,t)$
$\displaystyle=[1-\rho(x,t)]\int\rho_{A}(y,t)w_{\Gamma_{A}}([\phi_{A,B}],y,x,t)\,\mathrm{d}y$
(42)
$\displaystyle-\rho_{A}(x,t)\int[1-\rho(y,t)]w_{\Gamma_{A}}([\phi_{A,B}],x,y,t)\,\mathrm{d}y,$

and, by symmetry, a similar equation for $B$. The transition rates depend on
the utility function of each agent type and are a priori agent specific.
Denoting $u_{Z}(x)\equiv u_{Z}(x,[\phi_{A,B}])$ (with $Z=A$ or $B$), we set

$\displaystyle
w_{\Gamma_{Z}}([\phi_{A,B}],y,x,t)=\omega_{Z}f_{\Gamma_{Z}}[u_{Z}(x)-u_{Z}(y)],$
(43)

where $\omega_{Z}$ and $\Gamma_{Z}$ can be agent dependent.

Figure 7: Snapshots of the system for two frustrated interaction parameter
choices. (a) Stationary demixing in a region where LSA presents complex
eigenvalues. The agent $A$ phase still contains some $B$ agents. Parameters:
$c_{1}=-2$, $c_{2}=1$, $\sigma=3$, $\bar{\rho}_{A}=0.2$, $\bar{\rho}_{B}=0.5$,
$\Gamma=10$. (b) Chaotic propagation of polarized blobs in a region where LSA
presents pure real eigenvalues (null imaginary part). Parameters: $c_{1}=-2$,
$c_{2}=0.5$, $\sigma=7$, $\bar{\rho}_{A}=0.6$, $\bar{\rho}_{B}=0.2$,
$\Gamma=100$. For both (a) and (b), $L_{x}=L_{y}=300$. Movies are available
online [76].

In App. H, we perform the linear stability analysis of the homogeneous state.
As expected, in the frustrated two-population system, unstable modes can
display temporal oscillations. However, these oscillations may stop when
nonlinear terms become relevant, and the system may end up in a stationary
phase separation (similar to classical demixing in equilibrium systems), as
displayed in Fig. 7(a). Reciprocally, non-oscillating growing modes at the
linear level may give rise to propagating structures and waves when
nonlinearities become important, as shown in Fig. 7(b) (see Supplementary
Material in [76]). In our system, and at odds with recent work [74, 75], the
oscillatory nature of the non-homogeneous steady state cannot be predicted
from a simple linear stability analysis about the homogeneous solution.

A thorough analysis of the emerging behaviors in the multi-population system
would require more work, beyond the scope of the present paper. Still, it is
remarkable that, here as well, the linear stability analysis in the case of
local jumps yields exactly the same instability conditions as the nonlocal
dynamics ones (see results of Appendices H and I). As a consequence, we expect
that some results of the recent works [74, 75] should be relevant, to some
extent, to describe our multi-population system.

###  V.2 Housing market

A common and reasonable criticism of the kind of model developed here is that,
while the perceived density may be a significant factor in the decision making
process of agents, the price of a house should also necessarily be taken into
account. Indeed, in classical economics, the market is usually considered to
be the mechanism through which the optimal allocation of goods and assets
occurs (despite some contradicting empirical evidence e.g. [77]). As a result,
one could rightfully argue that a housing market is necessary to ensure that
agents eventually reach a steady state where their utility is maximal, at odds
with what we have observed in the condensed phase.

Incorporating pricing in the model is not trivial, however, and there are a
number of ways in which this could be done. A common approach in the modeling
of socioeconomic systems is to introduce an agent-dependent budget and to
constrain agents’ moves based on such budget, as done in [78] for example.
Realistically, this budget should then be heterogeneous among the agents (e.g.
Pareto distributed for instance). While relevant and interesting, this agent-
specific dependence as well as its formulation as a hard constraint would
require a different modeling approach, and is unlikely to be tractable
analytically. The alternative that we take here is to consider that if a given
move leads to excess costs for an agent, its utility would decrease. We may
then conveniently stay within the modeling framework of our model and assume
that such a price field $\psi$ is an increasing function of the smoothed
density field $\phi$, such that houses are more expensive in dense
neighborhoods– or in other words: prices are driven by demand only. In its
most general form, we propose the price-adjusted utility

$\bar{u}(\phi)=u(\phi)-u_{\mathrm{p}}[\psi(\phi)],$ (44)

where $u_{\mathrm{p}}$ is the price penalty, assumed to be an increasing
function of the price, and therefore, of the density $\phi$. This penalty is
then, by construction, expected to drive the system away from condensation.

For concreteness, one can consider the example of a linear penalty term on the
utility function , and that the price grows linearly with the local (smoothed)
density, such that:

$\bar{u}(\phi)=u(\phi)-\gamma\phi,$ (45)

with $\gamma>0$. Interestingly, introducing such a coupling boils down to
introducing a pair-wise repulsive interaction between agents. The condition to
observe condensation is then shifted by $\gamma$, and reads explicitly

$2\rho_{0}(1-\rho_{0})\frac{f^{\prime}_{\Gamma}(0)}{f_{\Gamma}(0)}(u^{\prime}(\rho_{0})-\gamma)>1.$
(46)

Clearly, condensation cannot occur anymore if the price penalty on the utility
is too important.

This example therefore illustrates that it is indeed possible for an
appropriate housing market to destroy the condensation observed in our
continuous model. However, it is important to note that the outcome (without
the price penalty) for agents is not necessarily improved by this brutal
homogenization through the price field. Besides, it can be argued that the
effect of price should not be as trivial as a utility penalty. In fact, other
models have used price as a proxy for social status, in which case agents are
actually more attracted by the most expensive site they can afford [78]. In
App. J we also consider the case of a more complex price field dependence.
Ongoing research is devoted to performing a comprehensive study of housing
price correlations in empirical data and the careful construction of an
adequate pricing dynamics [79].

##  VI Discussion

Let us summarize what we have achieved in this paper. We have introduced a
neighborhood-less extension of the Sakoda-Schelling model for the occupation
of a lattice representing a city. In this version of the model, the agents
attempt to maximize a utility that is a function of their perceived local
density, and are most satisfied when they are located in a site where such
density is of an intermediate value, i.e. neither empty nor too crowded.
Having that agents only consider their own site dependent utility, and that
their utility is nonlinear, drives the system out of equilibrium. As a result,
the system can no longer be studied by constructing a free energy directly
from an aggregate system-wide utility function, as was done by Grauwin et al.
[15] for agents inhabiting predefined neighborhoods or blocks in which the
utility is identical for all. Using numerical simulations as well as a mean-
field description of the nonequilibrium dynamics, we have established that the
apparent disparity between micromotives and macrobehaviours initially observed
by Schelling [4] is robust to the absence of neighborhoods and to the out-of-
equilibrium nature of our extension. Interestingly, we find that the
transition between the fully homogeneous state and the phase-separated one
likely belongs to the 2D Ising universality class, a debated topic in the
active matter literature regarding the very similar Motility Induced Phase
Separation (MIPS) phenomenon. Taking advantage of the similarity between our
problem and the Active Model B (describing MIPS), we predict the local density
in the bulk of the concentrated phase, confirming that the majority of agents
find themselves in a sub-optimal situations with a perceived density exceeding
the ideal value.

While seemingly technical, the fact that the original observations of
Schelling is robust to out-of-equilibrium dynamics actually carries far
reaching consequences, in our opinion. Indeed, as discussed in Sec. II.2,
equilibrium descriptions of socioeconomic problems require the decision rule
$f_{\Gamma}(x)$ to be the “logit” function. This very specific choice is a
common source of criticism, as any results are then a priori uncertain to hold
for other decision rules. Here, on the other hand, our out-of-equilibrium
description presents no such restriction, as all calculations have been
written as generally as possible and, interestingly, turn out to only depend
on $f_{\Gamma}(0)$ and $f^{\prime}_{\Gamma}(0)$. While the numerical results
presented here are those using the classical choice of the logit function for
lack of a more plausible decision rule, one could readily adapt the outcomes
following behavioral evidence that another function is more appropriate, and
we expect the entire phenomenology of our model to hold for any other sigmoid
function. More generally, we believe that this robustness to other decision
rules holds for a large number of socioeconomic models that have been
described using the methods of statistical physics [22, 80, 81, 82, 83]. Of
course, subtleties around the dynamics, such as the relaxation time towards
the steady state or the coarsening dynamics discussed here, will inherently be
affected by the specific choice that governs transition rates. This being
said, we have observed a remarkable similarity in the local and non-local
versions of our model for which the dynamics are yet qualitatively very
different. It is therefore possible that there also exists some degree of
universality in the dynamical behavior of different decision rules, at least
at the mean-field level.

Going back to the Sakoda-Schelling model, we have also introduced some
generalizations that we believe are natural and relevant. First, the
introduction of different sub-populations is interesting, as the exhibited
dynamical patterns are impossible to observe in an equilibrium version of the
model. Second, we have seen that introducing a linear dependence of the price
on the density has the effect of delaying the transition, eventually killing
it off completely if the price penalty in the utility function is strong
enough. As previously stated, however, this mechanism remains very simple. In
order to determine more plausible effects of a housing market, a thorough
analysis of real estate transactions appears to be a promising avenue, in
particular given the increasing availability of open datasets in this area in
major European cities. An extensive study of French data is currently
underway, hopefully allowing us to couple this continuous Sakoda-Schelling
model with a plausible housing market model in the near future [79]. Note
finally that the recent preprint [70] revisits the problem of urban
segregation in the United States and proposes a two-population model very
similar to ours. Studying the validity of hydrodynamic descriptions of the
two-population problem using the census data brought forth in this paper
appears to be an important perspective. Such a comparison could notably be
necessary to highlight the role of ingredients not taken into account in
existing models – such as public policy and economic inequalities – in the
emergence of strong geographic disparities.

## Acknowledgements

We warmly thank Jean-Philippe Bouchaud for his numerous insights on this
study, as well as Claire Alais and Noé Beserman who participated in the early
stages of this project. We also thank Eric Vanden-Eijnden for fruitful
discussions, as well as Eric Bertin for useful comments on the manuscript. R.
Z. also thanks Jérémy O’Byrne for precious discussions along the years on
thermodynamic mappings. J. G.-B. finally thanks Samy Lakhal for fruitful
discussions on the linear utility case. This research was conducted within the
Econophysics & Complex Systems Research Chair, under the aegis of the
Fondation du Risque, the Fondation de l’École polytechnique, the École
polytechnique and Capital Fund Management.

##  Appendix A Coarsening exponent

We start from a homogeneous system of size $L_{x}\times L_{y}$ with
$L_{x}=L_{y}$, and we quench it below the critical temperature in the spinodal
region. The system undergoes a spinodal decomposition where dense domains
coarsen until forming one single large cluster. The typical size of the
domains, denoted $L_{d}$ grows with time as $\sim t^{1/z}$, where the growth
exponent $z$ indicates the physics at play. To measure the typical domain
size, we compute first the structure factor, given by

$\displaystyle
S(\bm{k},t)=\Big{|}\sum_{\bm{r}}e^{-i\bm{k}\cdot\bm{r}}\phi(\bm{r},t)\Big{|}^{2}.$
(47)

Using isotropy of the system, we average the structure factor over given
shells $q=(k_{x}^{2}+k_{y}^{2})^{\frac{1}{2}}$ and we obtain the radial
structure factor $s(q,t)=\int_{[0,2\pi]}S(q,\theta,t)\mathrm{d}\theta$. The
typical domain size is given by

$\displaystyle
L_{d}(t)=2\pi\dfrac{\int_{k_{1}}^{\Lambda}s(q,t)\mathrm{d}q}{\int_{k_{1}}^{\Lambda}q\,s(q,t)\mathrm{d}q},$
(48)

with $\Lambda$ the ultraviolet cutoff and $k_{1}=2\pi/L_{x}$ the infrared
cutoff. On our finite grid, the integral takes the form of a discrete sum, the
wavenumber $q$ ranges from $2\pi/N_{x}$ to $2\pi(N_{x}-1)/N_{x}$, and the
increment $\mathrm{d}q$ is replaced by $2\pi/N_{x}$.

Figure 8: Snapshot of a Monte Carlo simulation. Green sites are occupied,
black sites are empty. We draw the boxes of size $\ell=L_{x}/6$ that are used
to measure the liquid and the gas densities. Here, the system size is
$L_{x}=200$, $L_{y}=66$. Parameters: $\alpha=3/2$, $\rho_{0}=0.35$, $T=0.05$.
Figure 9: Binder cumulant, order parameter and compressibility close to the
critical point $(\rho_{c},T_{c})=(0.271,0.0620)$ computed for $\alpha=3/2$ and
$\sigma=1$, as a function of the temperature $T$.

##  Appendix B Monte Carlo simulations for 2nd-order phase transitions

The Monte Carlo simulation setup that we used is shown in Fig. 8. As
introduced in recent works on MIPS, we study the phase transition using four
boxes of size $\ell\times\ell$ located in the bulks of the dense and “gas”
phases. The initial condition for the simulations is a fully separated state,
where a slab of density 1 coexists with a slab of density 0. We measure the
system decorrelation time $\tau_{d}$ and we start recording data after $\sim
1.5\tau_{d}$. Each simulation is run for a time $\geq 5\tau_{d}$, and each
symbol in Fig. 9 aggregates the data of $80$ independent simulations.

The collapse of the different observable with the 2D Ising critical exponents
is displayed in Fig. 9.

##  Appendix C Lyapunov function for non-local moves

We show below that, in the mean-field limit, with the logit decision rule and
a linear utility function, the hydrodynamic evolution has a Lyapunov function.
Starting from the pairwise Hamiltonian

$\mathcal{H}=-\frac{\nu}{2}\sum_{\bm{r},\bm{r}^{\prime}}n(\bm{r})G_{\sigma}(\bm{r}-\bm{r}^{\prime})n(\bm{r}^{\prime}),$
(49)

which can be shown to satisfy Eq. (5) when $u(\phi)=\nu\phi$, we take the
continuous-space limit and we account for the entropic contribution to find
the free energy functional

$\displaystyle\mathcal{F}[\rho]=-\frac{\nu}{2}\int\mathrm{d}x\,\mathrm{d}y\,\rho(x)G_{\sigma}(x-y)\rho(y)+TS[\rho],$
(50)

with the entropy $S[\rho]=\int[\rho\log(\rho)+(1-\rho)\log(1-\rho)]$. Before
addressing the dynamics of non-local moves, let us address the one of local
moves. It turns out that when moves are local, the use of a linear utility and
the logit rule implies that the mean-field dynamics is of gradient type,
analogous to a Wasserstein gradient flow but with $M[\rho]=\rho(1-\rho)$
instead of $M[\rho]=\rho$ [54], with $\mathcal{F}$ playing the role of the
free energy potential. Naturally, the stochastic dynamics is in detailed-
balance with the Gibbs-Boltzmann measure $e^{-\Gamma\mathcal{F}}$. When the
moves are no longer local, the gradient structure at the mean-field level is
difficult to unveil, even though the stochastic dynamics remains in detailed-
balance. The functional $\mathcal{F}$ remains a Lyapunov function of the
dynamics, i.e. it is non-increasing as the dynamics evolves. This is what we
show explicitly in the next paragraph.

For simplicity, we define the chemical potential

$\displaystyle\mu(x)\equiv\frac{\delta\mathcal{F}}{\delta\rho(x)}=-\nu\phi(x)+T\log\Big{(}\frac{\rho(x)}{1-\rho(x)}\Big{)},$
(51)

such that inverting the equation yields

$\displaystyle\rho(x)=\frac{e^{\frac{\Gamma}{2}(\mu(x)+\nu\phi(x))}}{D(x)},$
(52)

with $D(x)\equiv 2\cosh[\frac{\Gamma}{2}(\mu(x)+\nu\phi(x))]$. Then, computing
the total time derivative on the functional yields (we omit the time
dependence of the fields to alleviate notations):

$\displaystyle\begin{split}\frac{d\mathcal{F}}{dt}=&\int_{x}\frac{\delta\mathcal{F}}{\delta\rho(x)}\partial_{t}\rho(x,t)\mathrm{d}x\\\
=&\int_{x}\mu(x)(1-\rho(x,t))\int_{y}\rho(y,t)w_{\Gamma}([\phi],y,x,t)\mathrm{d}y\mathrm{d}x\\\
&-\int_{x}\mu(x)\rho(x,t)\int_{y}(1-\rho(y,t))w_{\Gamma}([\phi],x,y,t)\mathrm{d}y\mathrm{d}x\\\
=&-\iint\mathrm{d}x\mathrm{d}y\,\frac{Z(x,y)}{4D(x)D(y)\cosh[\frac{\Gamma}{2}(\phi(x)-\phi(y))]},\end{split}$
(53)

where we have replaced $w_{\Gamma}$ by the logit decision rate [see Eq. (16)]
and we have symmetrized and simplified the second line to obtain the third
line, and where we define
$Z(x,y)\equiv[\mu(y)-\mu(x)](e^{\frac{\Gamma}{2}(\mu(y)-\mu(x))}-e^{\frac{\Gamma}{2}(\mu(x)-\mu(y))})\geq
0$, $\forall x,y$. We thus conclude that $\frac{d\mathcal{F}}{dt}\leq 0$, i.e.
that $\mathcal{F}$ is a Lyapunov function of the hydrodynamic evolution when
the utility is linear. Now, more than that, the function is always strictly
decreasing unless it starts from a fixed point, which forbids limit cycles.
Indeed, one notices that the integrand in Eq. (53) is always positive, unless
$\mu(x,t)=\mu(y,t)$, for all $x,y$. Setting $C(t)=\mu(x,t)$ on the manifold
where $\mathcal{F}$ is constant, we have

$\displaystyle 1-\rho(x,t)=\rho(x,t)e^{\Gamma C(t)}e^{\Gamma\nu\phi(x,t)}.$
(54)

If we inject this relation in the non-local mean-field evolution Eq. (15),
then we obtain $\partial_{t}\rho(x,t)=0$, indicating that $\rho(x,t)$ must be
a stationary fixed point when the Lyapunov function is constant.

##  Appendix D Local mean field description and LSA

In this section, we consider a modified dynamics where agents are allowed to
relocate on neighboring site only. For simplicity, we also consider that the
system is one dimensional. It is thus possible to perform a Taylor expansion
of the different fields assuming that all fields are smooth in the mean-field
limit. The jump probability between two neighboring sites becomes

$\displaystyle
f_{\Gamma}[u(x+a)-u(x)]=f_{\Gamma}\left(a\partial_{x}u+\frac{a^{2}}{2}\partial_{x}^{2}u\right),$
(55)

where $a$ is the lattice size, and $u$ is the utility on position $x$. The
evolution of the density (for non-overlapping agents) is thus given by

$\displaystyle\begin{split}\partial_{t}\rho=&\rho(x+a)[1-\rho(x)]f_{\Gamma}(-a\partial_{x}u-\frac{a^{2}}{2}\partial_{x}^{2}u)\\\
&+\rho(x-a)[1-\rho(x)]f_{\Gamma}(a\partial_{x}u-\frac{a^{2}}{2}\partial_{x}^{2}u)\\\
&-\rho(x)[1-\rho(x+a)]f_{\Gamma}(a\partial_{x}u+\frac{a^{2}}{2}\partial_{x}^{2}u)\\\
&-\rho(x)[1-\rho(x-a)]f_{\Gamma}(-a\partial_{x}u+\frac{a^{2}}{2}\partial_{x}^{2}u).\end{split}$
(56)

After Taylor expansion up to $O(a^{2})$ and time rescaling, it turns out that
the evolution equation simplifies into

$\displaystyle\partial_{t}\rho$
$\displaystyle=f_{\Gamma}(0)\partial_{x}^{2}\rho-2f_{\Gamma}^{\prime}(0)\partial_{x}[\rho(1-\rho)\partial_{x}u].$
(57)

Then, expanding around an homogeneous state, we write
$\rho=\rho_{0}+\rho_{1}(x,t)$, $\phi=\rho_{0}+\phi_{1}(x,t)$, and we obtain to
leading order in the perturbation:

$\partial_{t}\rho_{1}=f_{\Gamma}(0)\partial_{x}^{2}\rho_{1}-2f_{\Gamma}^{\prime}(0)\rho_{0}(1-\rho_{0})u^{\prime}(\rho_{0})\partial_{x}^{2}\phi_{1}.$
(58)

In Fourier space the evolution of the mode $k$ is given by
$\partial_{t}\hat{\rho}_{1}=\Lambda(k)\hat{\rho}_{1}$, with

$\Lambda(k)=-k^{2}f_{\Gamma}(0)\left(1-2\frac{f^{\prime}_{\Gamma}(0)}{f_{\Gamma}(0)}\rho_{0}(1-\rho_{0})u^{\prime}(\rho_{0})\hat{G}_{\sigma}(k)\right).$
(59)

From this, we deduce that the homogeneous system is unstable if there exists a
mode $k^{\star}$ such that

$\displaystyle
1<2\frac{f^{\prime}_{\Gamma}(0)}{f_{\Gamma}(0)}\rho_{0}(1-\rho_{0})u^{\prime}(\rho_{0})\hat{G}_{\sigma}(k^{\star}).$
(60)

This criterion is exactly the same as the one found for the non-local move
dynamics.

##  Appendix E Local versus non-local PDEs

To illustrate the effectiveness of our local-move approximation in the
description of the steady state of the system, we have solved numerically both
the local and the non-local mean-field PDEs for the same parameters. The
resulting density profiles, displayed in Fig. 10, appear to be strictly
identical up to numerical errors.

Figure 10: 1D steady-state profiles of the density $\rho$ computed by solving
the mean-field dynamics with local moves (solid line), or non-local moves
(dashed line). The two density profiles superimpose almost exactly,
independently of the various parameter values. It was checked that different
initial configurations lead to the same final state. a) Parameters: utility
exponent $\alpha=1$, $\rho_{0}=0.4$, $\sigma=12$, $L=300$, $\Gamma=9$. b)
$\alpha=3/2$, $\rho_{0}=0.3$, $\sigma=10$, $L=300$, $\Gamma=15$.

##  Appendix F Breakdown of the gradient expansion

To illustrate the possible breakdown of the gradient expansion even for
$\kappa(\rho)>0$, we consider a strictly increasing and continuous utility
with a $|\rho-\rho^{\star}|^{-1/2}$ divergence. As shown in the comparison
with the mean-field PDE solved numerically in Fig. 11, the generalized
thermodynamics fails at predicting the bulk densities in this pathological
case.

Figure 11: a) Monotonic utility
$u(\rho)=\operatorname*{sgn}(\rho-\rho^{\star})|\rho-\rho^{\star}|^{1/2}$
(solid line), and bijective change of variable
$R(\rho)=\operatorname*{sgn}(\rho-\rho^{\star})|\rho-\rho^{\star}|^{3/2}$
(dashed line) for $\rho^{\star}=1/2$. b) Density profile when phase separation
occurs, for $\Gamma=2.5$, $\sigma=10$, $L_{x}=300$. The plateaus of the liquid
and the gas phase do not match the plateaus predicted by the generalized
thermodynamic mapping because the gradient expansion is no longer valid close
to $\rho=0.5$.

##  Appendix G Details on the double-tangent and the Maxwell constructions

For completeness, we provide a summary of the two approaches that yields the
binodal densities, as they are presented in [65]. As stated the main text,
when a phase separation occurs in equilibrium, the density profile in the
stationary state is no longer evolving, i.e. the free energy, constrained by
the fact that the total mass of the field is fixed, has reached a minimum.
Since interfaces have a sub-extensive contribution in the thermodynamics
limit, one can work with the free energy density $f(\rho)$. The chemical
potential it thus $\mu(\rho)=\frac{\mathrm{d}f}{\mathrm{d}\rho}$. In the
steady state, one has equality of the chemical potential in the liquid and in
the gas, i.e

$\displaystyle\mu(\rho_{\ell})=\mu(\rho_{g})=\bar{\mu},$ (61)

where $\rho_{\ell}$ and $\rho_{g}$ denote the densities in the liquid and in
the gas, respectively. Since the interface does not move, one also has
equality of pressure across the interface, i.e

$\displaystyle P(\rho_{\ell})=P(\rho_{g})=\bar{P},$ (62)

with $P(\rho)=\rho\mu(\rho)-f(\rho)$. We are thus looking for the two
densities such that the tangent to the free energy density $f$ in
$\rho_{\ell}$ and $\rho_{g}$ is the same, with the slope given by

$\displaystyle\bar{\mu}=\frac{f(\rho_{\ell})-f(\rho_{g})}{\rho_{\ell}-\rho_{g}}.$
(63)

Equivalently, the coexisting densities can be obtained via the Maxwell equal-
area construction that imposes

$\displaystyle\int_{v_{\ell}}^{v_{g}}[P(v)-\bar{P}]\mathrm{d}v=0,$ (64)

where $v\equiv 1/\rho$ is the volume per particle, and
$v_{g/\ell}=1/\rho_{g/\ell}$.

Figure 12: (a) Double-tangent construction (in red) on the function $\Phi(R)$
and (b) Maxwell equal-area construction on $h_{0}(w)$, with $h_{0}(w)=\bar{h}$
in red. Shaded areas in (b) are equal. Here we have plotted these functions
for $\alpha=3/2$, $\rho^{\star}=1/2$, $f_{\Gamma}(0)=1/2$,
$f^{\prime}_{\Gamma}(0)=\Gamma/4$, and $\Gamma=10$. We find $\bar{h}=0.0212$,
$R_{\ell}=1.092$ and $R_{g}=0.327$, (or $w_{\ell}=0.915$ and $w_{g}=3.058$),
translating into $\rho_{\ell}=0.508$ and $\rho_{g}=0.047$.

Here, even though we lie out of equilibrium, we have shown how to find a
function $\Phi(R)$ that plays a role similar to a free energy density by means
of the change of variable $\rho(R)$. This function is always convex for $T$
above the mean-field critical temperature $T_{c}^{\mathrm{MF}}$ and display a
non-convex region for $T<T_{c}^{\mathrm{MF}}$. The function
$g_{0}(R)=\frac{\mathrm{d}\Phi}{\mathrm{d}R}$ is now the chemical potential,
and the double construction on $\Phi$ imposes

$\displaystyle
g_{0}(R_{\ell})=g_{0}(R_{g})=\frac{\Phi(R_{\ell})-\Phi(R_{g})}{R_{\ell}-R_{g}}.$
(65)

The equal-area Maxwell construction involves the so-called generalized
pressure $h_{0}$:

$\displaystyle h_{0}(R)=R\frac{\mathrm{d}\phi(R)}{\mathrm{d}R}-\Phi(R),$ (66)

and setting $w=1/R$, we find $w_{\ell}$ and $w_{g}$ such that

$\displaystyle\int_{w_{\ell}}^{w_{g}}[h_{0}(w)-\bar{h}]\mathrm{d}w=0.$ (67)

In practice, the function $h_{0}$ can be obtained numerically and the volume
$w_{\ell}$ and $w_{g}$ can be obtained with a numerical solver using Eq. (67),
more easily than solving the double-tangent condition. We show both
constructions in Fig. 12. Note again that $R$ has no clear physical meaning
and is an intermediary variable of computations. As such, one can always
choose integration constants in Eq. (35) such that $R(\rho)\geq R(0)>0$ to
ensure that $w=1/R$ is correctly defined. In the main text though, we have
chosen $R(\rho)$ centered in $0$ for simplicity, and computed $g_{0}(R)$ and
$\Phi(R)$ accordingly. If we enforce below $R(0)=1-\rho^{\star}{}^{2-\alpha}$,
then taking $\alpha=3/2$, $\rho^{\star}=1/2$, $f_{\Gamma}(0)=1/2$ and
$f^{\prime}_{\Gamma}(0)=\Gamma/4$, we can obtain an analytic expression for
$\Phi(R)$

$\displaystyle\begin{split}\Phi(R)=&\frac{\sqrt{2}}{2}\left(2\Theta(R-1)-1\right)\left[\frac{R-1}{\sqrt{2}}\log\left(\frac{1+2(R-1)^{2}}{1-2(R-1)^{2}}\right)+\arctan\left((R-1)\sqrt{2}\right)-\operatorname{arctanh}\left((R-1)\sqrt{2}\right)\right]\\\
&+\frac{\Gamma}{8}|R-1|(R-1)^{3}.\end{split}$ (68)

$\displaystyle\begin{split}\Phi(R)=\sqrt{2}(2\Theta(R-1)-1)\left[\frac{R-1}{2\sqrt{2}}\log\left(\frac{2+(R-1)^{2}}{2-(R-1)^{2}}\right)+\arctan\left(\frac{R-1}{\sqrt{2}}\right)-\operatorname{arctanh}\left(\frac{R-1}{\sqrt{2}}\right)\right]+\frac{\Gamma}{64}|R-1|(R-1)^{3}.\end{split}$
(69)

$\displaystyle\begin{split}\Phi(R)=&(2\Theta(R)-1)\left[\frac{R}{2}\log\left(\frac{2+R^{2}}{2-R^{2}}\right)+\sqrt{2}\arctan\left(\frac{R}{\sqrt{2}}\right)-\sqrt{2}\operatorname{arctanh}\left(\frac{R}{\sqrt{2}}\right)\right]+\frac{1}{64}\Gamma|R|R^{3}\end{split}$
(70)

##  Appendix H Linear stability for two coupled populations

We consider the evolution of a perturbation of the homogeneous state in Eq.
(42) (and in its coupled analogue for the field $\rho_{B}$). Close to the
homogeneous state $\rho_{A}(x)\equiv\bar{\rho}_{A}$,
$\rho_{B}(x)\equiv\bar{\rho}_{B}$, with
$\rho_{0}=\bar{\rho}_{A}+\bar{\rho}_{B}$, we expand the fields
$\rho_{Z}(x,t)=\bar{\rho}_{Z}+\rho_{Z,1}(x,t)$, with $Z=A$ or $B$, and the
perturbation fields are denoted with index $1$. One also has
$\rho(x,t)=\rho_{0}+\rho_{1}(x,t)$, and
$\phi_{Z}(x,t)=\bar{\rho}_{Z}+\phi_{Z,1}(x,t)$. Keeping leading order terms in
Eq. (42) yields

$\displaystyle\partial_{t}\rho_{A,1}=\Omega\omega_{A}\left[-\bar{\rho}_{A}f_{\Gamma_{A}}(0)\rho_{1}(x,t)+2(1-\rho_{0})\bar{\rho}_{A}f^{\prime}_{\Gamma_{A}}(0)[\phi_{A,1}(x,t)\partial_{1}u_{A}+\phi_{B,1}(x,t)\partial_{2}u_{A}]-(1-\rho_{0})f_{\Gamma_{A}}(0)\rho_{A,1}(x,t)\right],$
(71)

where $\partial_{1}u_{A}$ is a shorthand notation for $\frac{\partial
u_{A}}{\partial\bar{\rho}_{A}}[\bar{\rho}_{A},\bar{\rho}_{B}]$. Taking the
logistic function $f_{\Gamma_{A}}(0)=\frac{1}{2}$,
$f^{\prime}_{\Gamma_{A}}(0)=\frac{\Gamma_{A}}{4}$, the linear evolution
simplifies into

$\displaystyle\partial_{t}\rho_{A,1}(x,t)=\frac{\Omega\omega_{A}}{2}\left[-\bar{\rho}_{A}\rho_{1}(x,t)+(1-\rho_{0})\bar{\rho}_{A}\Gamma_{A}[\phi_{A,1}(x,t)\partial_{1}u_{A}+\phi_{B,1}(x,t)\partial_{2}u_{A}]-(1-\rho_{0})\rho_{A,1}(x,t)\right].$
(72)

Similarly, we obtain for the evolution of $B$:

$\displaystyle\partial_{t}\rho_{B,1}(x,t)=\frac{\Omega\omega_{B}}{2}\left[-\bar{\rho}_{B}\rho_{1}(x,t)+(1-\rho_{0})\bar{\rho}_{B}\Gamma_{B}[\phi_{A,1}(x,t)\partial_{1}u_{B}+\phi_{B,1}(x,t)\partial_{2}u_{B}]-(1-\rho_{0})\rho_{B,1}(x,t)\right].$
(73)

Denoting $\hat{\rho}_{Z}(k,t)$ the Fourier transform of $\rho_{Z,1}(x,t)$, the
evolution equation can be cast in Fourier space into

$\displaystyle\partial_{t}\begin{pmatrix}\hat{\rho}_{A}(k,t)\\\
\hat{\rho}_{B}(k,t)\end{pmatrix}=L\begin{pmatrix}\hat{\rho}_{A}(k,t)\\\
\hat{\rho}_{B}(k,t)\end{pmatrix},$ (74)

with

$\displaystyle
L=\frac{\Omega}{2}\begin{pmatrix}\omega_{A}(\bar{\rho}_{B}-1+(1-\rho_{0})\bar{\rho}_{A}\Gamma_{A}\hat{G}_{\sigma}(k)\partial_{1}u_{A})&\omega_{A}(-\bar{\rho}_{A}+(1-\rho_{0})\bar{\rho}_{A}\Gamma_{A}\hat{G}_{\sigma}(k)\partial_{2}u_{A})\\\
\omega_{B}(-\bar{\rho}_{B}+(1-\rho_{0})\bar{\rho}_{B}\Gamma_{B}\hat{G}_{\sigma}(k)\partial_{1}u_{B})&\omega_{B}(\bar{\rho}_{A}-1+(1-\rho_{0})\bar{\rho}_{B}\Gamma_{B}\hat{G}_{\sigma}(k)\partial_{2}u_{B})\end{pmatrix}.$
(75)

For simplicity, we will consider that agents are equally rational
($\Gamma_{A}=\Gamma_{B}=\Gamma$) and that their moving rates are also
identical ($\omega_{A}=\omega_{B}=\omega$).

We are looking for conditions to observe dynamical patterns and/or static
phase separation. Notably, the homogeneous state is linearly unstable if one
eigenvalue of $L$ has a positive real part. It is important to stress that the
linear stability analysis is unable to predict the dynamic behavior when
nonlinear terms become relevant. Whether the eigenvalues display an imaginary
part or not _does not_ bring any information on the final dynamics of the
system. For the sake of completeness, we explicitate the criteria to have
eigenvalues with positive real part and zero imaginary part, referred to as
case (i), and eigenvalues with positive real part and nonzero imaginary part,
referred to as case (ii). We lie in case (i) if

$\displaystyle\begin{cases}\operatorname*{tr}L>0\\\
(\operatorname*{tr}L)^{2}-4\det L>0,\end{cases}\text{ or
}\begin{cases}\operatorname*{tr}L<0\\\ \det L<0.\end{cases}$ (76)

Case (ii) is obtained if

$\displaystyle\begin{cases}\operatorname*{tr}L>0\\\
(\operatorname*{tr}L)^{2}-4\det L<0.\end{cases}$ (77)

The criterion $\operatorname*{tr}L>0$ notably simplifies into

$\displaystyle\bar{\rho}_{A}\partial_{1}u_{A}+\bar{\rho}_{B}\partial_{2}u_{B}>\frac{1}{\Gamma\hat{G}_{\sigma}(k)}\left(\frac{2-\rho_{0}}{1-\rho_{0}}\right).$
(78)

In the main text we have come up with utility functions that lead to
eigenvalues with positive real parts and non zero imaginary parts, thus
suggesting chasing instability. In some cases, oscillations were observed
close to the homogeneous state but they eventually vanished at late times.
Whether or not the chasing instability or oscillations are sustained cannot be
predicted from the simple linear stability analysis but would require to
perform a weakly non-linear analysis which is beyond the scope of this present
paper.

##  Appendix I LSA for two populations with local moves

We start from the local jump approximation of the mean-field equation for the
coupled fields. We find that the dynamics can be cast into

$\partial_{t}\rho_{A}=\partial_{x}[\rho_{A}(1-\rho_{A}-\rho_{B})\partial_{x}\mu([\rho_{A,B}],x)],$
(79)

with $\mu=\mu_{\mathrm{ent.}}+\mu_{\mathrm{util.}}$,

$\mu_{\mathrm{ent.}}=w_{\Gamma_{A}}(0)\log\left(\frac{\rho_{A}}{1-\rho_{A}-\rho_{B}}\right),$
(80) $\mu_{\mathrm{util.}}=-2w_{\Gamma_{A}}^{\prime}(0)u_{A}([\rho],x),$ (81)

and likewise for $\rho_{B}$. One can look into the stability of an homogeneous
state with densities $\bar{\rho}_{A}$ and $\bar{\rho}_{B}$, expanding around
this state with a utility $u(\phi_{A},\phi_{B})$ for agents $A$ and
$v(\phi_{A},\phi_{B})$ for agents $B$. For convenience, we will take
$w_{\Gamma_{A}}(0)=w_{\Gamma_{B}}(0)=\omega f_{\Gamma}(0)=\omega/2$ and
$w_{\Gamma_{A}}^{\prime}(0)=w_{\Gamma_{B}}^{\prime}(0)=\omega
f_{\Gamma}^{\prime}(0)=\omega\Gamma/4$. Expanding around the homogeneous state
$(\bar{\rho}_{A}$ , $\bar{\rho}_{B})$ leads to

$\begin{cases}\partial_{t}\rho_{A,1}=\frac{\omega}{2}\left[(1-\bar{\rho}_{B})\partial_{x}^{2}\rho_{A,1}+\bar{\rho}_{A}\partial_{x}^{2}\rho_{B,1}-\partial_{x}(\Gamma\bar{\rho}_{A}(1-\rho_{0})\partial_{x}\phi_{A,1}\partial_{1}u+\partial_{x}\phi_{B,1}\partial_{2}u)\right]\\\
\partial_{t}\rho_{B,1}=\frac{\omega}{2}\left[(1-\bar{\rho}_{A})\partial_{x}^{2}\rho_{B,1}+\bar{\rho}_{B}\partial_{x}^{2}\rho_{A,1}-\partial_{x}(\Gamma\bar{\rho}_{B}(1-\rho_{0})\partial_{x}\phi_{A,1}\partial_{1}v+\partial_{x}\phi_{B,1}\partial_{2}v)\right],\end{cases}$
(82)

Hence, in Fourier space, the linear system can be cast into

$\displaystyle\partial_{t}\begin{pmatrix}\hat{\rho}_{A}(k,t)\\\
\hat{\rho}_{B}(k,t)\end{pmatrix}=K\begin{pmatrix}\hat{\rho}_{A}(k,t)\\\
\hat{\rho}_{B}(k,t)\end{pmatrix},$ (83)

with

$\displaystyle K=\frac{\omega
k^{2}}{2}\begin{pmatrix}\bar{\rho}_{B}-1+\Gamma\bar{\rho}_{A}(1-\rho_{0})\hat{G}_{\sigma}(k)\partial_{1}u&-\bar{\rho}_{A}+\Gamma\bar{\rho}_{A}(1-\rho_{0})\hat{G}_{\sigma}(k)\partial_{2}u\\\
-\bar{\rho}_{B}+(1-\rho_{0})\bar{\rho}_{B}\Gamma\hat{G}_{\sigma}(k)\partial_{1}v&\bar{\rho}_{A}-1+(1-\rho_{0})\bar{\rho}_{B}\Gamma\hat{G}_{\sigma}(k)\partial_{2}v\end{pmatrix}.$
(84)

It is interesting to note that the evolution matrix $K$ is directly
proportional to $L$ and, as a consequence, the stability criterion of the
homogeneous state with local moves is exactly the same as the one found for
non-local moves.

##  Appendix J Introducing a non-linear price field

Figure 13: a) Non-monotonic price field $\psi(\phi)$ as a function of the
locally perceived density $\phi$, given by Eq. (85). Parameters:
$\rho_{\mathrm{p}}^{\star}=0.25$ and ${\alpha_{\mathrm{p}}}=2$. b) Spinodal
curves for parameters $\lambda=5,{\alpha_{\mathrm{p}}}=2,\alpha=3$ and
densities [$\rho^{\star}=0.2,\rho_{\mathrm{p}}^{\star}=0.4$] (solid line),
[$\rho^{\star}=0.45,\rho_{\mathrm{p}}^{\star}=0.2$] (dotted line) and
[$\rho^{\star}=0.5,\rho_{\mathrm{p}}^{\star}=0.1$] (dashed line).

For completeness, we can also consider a non-monotonic price field (see Fig.
13(a))

$\psi=\rho_{\mathrm{p}}^{\star}-\lvert\phi-\rho_{\mathrm{p}}^{\star}\rvert^{{\alpha_{\mathrm{p}}}}.$
(85)

The intuition behind this relation is that prices can be lower in overcrowded
neighborhoods as well as in empty neighborhoods, and are maximized for a given
density $\rho_{\mathrm{p}}^{\star}$. Assuming again that the price-adjusted
utility has the form $\bar{u}(\phi)=u(\phi)-u_{\mathrm{p}}[\psi(\phi)]$, with
$u_{\mathrm{p}}[\psi]=\lambda\psi$ and $\lambda>0$, the total utility for the
agents is then given by

$\bar{u}(\phi)=-\lvert\phi-\rho^{\star}\rvert^{\alpha}+\lambda\lvert\phi-\rho_{\mathrm{p}}^{\star}\rvert^{\alpha_{\mathrm{p}}}+cst,$
(86)

We can inject this expression into the linear-stability condition, see Eq.
(21), to pinpoint the condensation. In Fig. 13(b) we take $\alpha=3$,
$\lambda=5,{\alpha_{\mathrm{p}}}=2$, and we interestingly observe that some
spinodal curves display several re-entrance points.

## References

  * Hotblling [1929] H. Hotblling, Stability in competition, The economic journal 39, 41 (1929).
  * Pigou [2017] A. Pigou,  _The economics of welfare_ (Routledge, 2017).
  * Braess [1968] D. Braess, Über ein paradoxon aus der verkehrsplanung, Unternehmensforschung 12, 258 (1968).
  * Schelling [1971] T. C. Schelling, Dynamic models of segregation, Journal of mathematical sociology 1, 143 (1971).
  * Hegselmann [2017] R. Hegselmann, Thomas C. Schelling and James M. Sakoda: The Intellectual, Technical, and Social History of a Model, Journal of Artificial Societies and Social Simulation 20, 15 (2017).
  * Sakoda [1971] J. M. Sakoda, The checkerboard model of social interaction, The Journal of Mathematical Sociology 1, 119 (1971).
  * Boustan [2013] L. P. Boustan,  _Racial residential segregation in American cities_ , Tech. Rep. (National Bureau of Economic Research, 2013).
  * Trounstine [2018] J. Trounstine,  _Segregation by design: Local politics and inequality in American cities_ (Cambridge University Press, 2018).
  * Schelling [2006] T. C. Schelling,  _Micromotives and macrobehavior_ (WW Norton & Company, 2006).
  * Vinković and Kirman [2006] D. Vinković and A. Kirman, A physical analogue of the Schelling model, Proceedings of the National Academy of Sciences 103, 19261 (2006), publisher: Proceedings of the National Academy of Sciences.
  * Dall’Asta  _et al._ [2008] L. Dall’Asta, C. Castellano, and M. Marsili, Statistical physics of the Schelling model of segregation, Journal of Statistical Mechanics: Theory and Experiment 2008, L07002 (2008).
  * Gauvin  _et al._ [2009] L. Gauvin, J. Vannimenus, and J.-P. Nadal, Phase diagram of a Schelling segregation model, The European Physical Journal B 70, 293 (2009).
  * Rogers and McKane [2011] T. Rogers and A. J. McKane, A unified framework for Schelling’s model of segregation, Journal of Statistical Mechanics: Theory and Experiment 2011, P07006 (2011).
  * Gauvin  _et al._ [2010] L. Gauvin, J.-P. Nadal, and J. Vannimenus, Schelling segregation in an open city: A kinetically constrained Blume-Emery-Griffiths spin-1 system, Physical Review E 81, 066120 (2010).
  * Grauwin  _et al._ [2009] S. Grauwin, E. Bertin, R. Lemoy, and P. Jensen, Competition between collective and individual dynamics, Proceedings of the National Academy of Sciences 106, 20622 (2009).
  * Jensen  _et al._ [2018] P. Jensen, T. Matreux, J. Cambe, H. Larralde, and E. Bertin, Giant Catalytic Effect of Altruists in Schelling’s Segregation Model, Physical Review Letters 120, 208301 (2018).
  * Barmpalias  _et al._ [2015] G. Barmpalias, R. Elwes, and A. Lewis-Pye, Tipping points in 1-dimensional Schelling models with switching agents, Journal of Statistical Physics 158, 806 (2015).
  * Ortega  _et al._ [2021a] D. Ortega, J. Rodríguez-Laguna, and E. Korutcheva, A Schelling model with a variable threshold in a closed city segregation model. analysis of the universality classes, Physica A: Statistical Mechanics and its Applications 574, 126010 (2021a).
  * Ortega  _et al._ [2021b] D. Ortega, J. Rodríguez-Laguna, and E. Korutcheva, Avalanches in an extended Schelling model: An explanation of urban gentrification, Physica A: Statistical Mechanics and its Applications 573, 125943 (2021b).
  * Abella  _et al._ [2022] D. Abella, M. San Miguel, and J. J. Ramasco, Aging effects in Schelling segregation model, Scientific Reports 12, 19376 (2022).
  * Bouchaud [2013] J.-P. Bouchaud, Crises and Collective Socio-Economic Phenomena: Simple Models and Challenges, Journal of Statistical Physics 151, 567 (2013).
  * Garnier-Brun  _et al._ [2023] J. Garnier-Brun, J.-P. Bouchaud, and M. Benzaquen, Bounded rationality and animal spirits: A fluctuation-response approach to Slutsky matrices, Journal of Physics: Complexity 4 (2023).
  * Wittkowski  _et al._ [2014] R. Wittkowski, A. Tiribocchi, J. Stenhammar, R. J. Allen, D. Marenduzzo, and M. E. Cates, Scalar $\varphi^{4}$ field theory for active-particle phase separation, Nature Communications 5, 4351 (2014).
  * Cahn and Hilliard [1958] J. W. Cahn and J. E. Hilliard, Free energy of a nonuniform system. i. interfacial free energy, Journal of Chemical Physics 28, 258 (1958).
  * Luce [1977] R. D. Luce, The choice axiom after twenty years, Journal of mathematical psychology 15, 215 (1977).
  * Anderson  _et al._ [1992] S. P. Anderson, A. De Palma, and J.-F. Thisse,  _Discrete choice theory of product differentiation_ (MIT press, 1992).
  * Luce [1959] R. D. Luce,  _Individual choice behavior_ (Wesley, 1959).
  * Nadal  _et al._ [1998] J.-P. Nadal, G. Weisbuch, O. Chenevez, and A. Kirman, A formal approach to market organization: choice functions, mean field approximation and maximum entropy principle, Advances in Self-Organization and Evolutionary Economics , 149 (1998).
  * Marsili [1999] M. Marsili, On the multinomial logit model, Physica A: Statistical Mechanics and its Applications 269, 9 (1999).
  * Brock and Durlauf [2001] W. A. Brock and S. N. Durlauf, Discrete choice with social interactions, The Review of Economic Studies 68, 235 (2001).
  * Grauwin  _et al._ [2012] S. Grauwin, F. Goffette-Nagot, and P. Jensen, Dynamic models of residential segregation: An analytical solution, Journal of Public Economics 96, 124 (2012).
  * Kolmogoroff [1936] A. Kolmogoroff, Zur Theorie der Markoffschen Ketten, Mathematische Annalen 112, 155 (1936).
  * Tamayo and Klein [1989] P. Tamayo and W. Klein, Critical dynamics and global conservation laws, Physical Review Letters 63, 2757 (1989).
  * Bray [1991] A. J. Bray, Comment on “Critical dynamics and global conservation laws”, Physical Review Letters 66, 2048 (1991).
  * Rutenberg [1996] A. D. Rutenberg, Nonequilibrium phase ordering with a global conservation law, Physical Review E 54, 972 (1996).
  * Tartaglia  _et al._ [2018] A. Tartaglia, L. F. Cugliandolo, and M. Picco, Coarsening and percolation in the kinetic 2 d Ising model with spin exchange updates and the voter model, Journal of Statistical Mechanics: Theory and Experiment 2018, 083202 (2018).
  * Siebert  _et al._ [2018] J. T. Siebert, F. Dittrich, F. Schmid, K. Binder, T. Speck, and P. Virnau, Critical behavior of active brownian particles, Physical Review E 98, 030601 (2018).
  * Maggi  _et al._ [2021] C. Maggi, M. Paoluzzi, A. Crisanti, E. Zaccarelli, and N. Gnan, Universality class of the motility-induced critical point in large scale off-lattice simulations of active particles, Soft Matter 17, 3807 (2021).
  * Dittrich  _et al._ [2021] F. Dittrich, T. Speck, and P. Virnau, Critical behavior in active lattice models of motility-induced phase separation, The European Physical Journal E 44, 53 (2021).
  * Rovere  _et al._ [1990] M. Rovere, D. W. Heermann, and K. Binder, The gas-liquid transition of the two-dimensional Lennard-Jones fluid, Journal of Physics: Condensed Matter 2, 7009 (1990).
  * Rovere  _et al._ [1993] M. Rovere, P. Nielaba, and K. Binder, Simulation studies of gas-liquid transitions in two dimensions via a subsystem-block-density distribution analysis, Zeitschrift für Physik B Condensed Matter 90, 215 (1993).
  * Binder [1997] K. Binder, Applications of monte carlo methods to statistical physics, Reports on Progress in Physics 60, 487 (1997).
  * Gnan and Maggi [2022] N. Gnan and C. Maggi, Critical behavior of quorum-sensing active particles, Soft Matter 18, 7654 (2022).
  * Nakano and Adachi [2023] H. Nakano and K. Adachi, Universal properties of repulsive self-propelled particles and attractive driven particles, arXiv preprint arXiv:2306.17517 (2023).
  * Caballero  _et al._ [2018] F. Caballero, C. Nardini, F. van Wijland, and M. E. Cates, Strong Coupling in Conserved Surface Roughening: A New Universality Class?, Physical Review Letters 121, 020601 (2018).
  * Speck [2022] T. Speck, Critical behavior of active Brownian particles: Connection to field theories, Physical Review E 105, 064601 (2022).
  * Curatolo  _et al._ [2020] A. I. Curatolo, N. Zhou, Y. Zhao, C. Liu, A. Daerr, J. Tailleur, and J. Huang, Cooperative pattern formation in multi-component bacterial systems through reciprocal motility regulation, Nature Physics 16, 1152 (2020).
  * Curatolo  _et al._ [2016] A. I. Curatolo, M. R. Evans, Y. Kafri, and J. Tailleur, Multilane driven diffusive systems, Journal of Physics A: Mathematical and Theoretical 49, 095601 (2016), publisher: IOP Publishing.
  * Allen and Cahn [1979] S. M. Allen and J. W. Cahn, A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta metallurgica 27, 1085 (1979).
  * Cahn and Hilliard [2004] J. W. Cahn and J. E. Hilliard, Free Energy of a Nonuniform System. I. Interfacial Free Energy, The Journal of Chemical Physics 28, 258 (2004).
  * Grafke  _et al._ [2017] T. Grafke, M. E. Cates, and E. Vanden-Eijnden, Spatiotemporal Self-Organization of Fluctuating Bacterial Colonies, Physical Review Letters 119, 188003 (2017).
  * O’Byrne and Tailleur [2020] J. O’Byrne and J. Tailleur, Lamellar to Micellar Phases and Beyond: When Tactic Active Systems Admit Free Energy Functionals, Physical Review Letters 125, 208003 (2020).
  * O’Byrne [2023] J. O’Byrne, Nonequilibrium currents in stochastic field theories: A geometric insight, Physical Review E 107, 054105 (2023).
  * Otto [2001] F. Otto, The Geometry of Dissipative Evolution Equations: The Porous Medium Equation, Communications in Partial Differential Equations 26, 101 (2001), publisher: Taylor & Francis _eprint: https://doi.org/10.1081/PDE-100002243.
  * Lefevre and Biroli [2007] A. Lefevre and G. Biroli, Dynamics of interacting particle systems: stochastic process and field theory, Journal of Statistical Mechanics: Theory and Experiment 2007, P07024 (2007).
  * Tailleur and Cates [2008] J. Tailleur and M. E. Cates, Statistical Mechanics of Interacting Run-and-Tumble Bacteria, Physical Review Letters 100, 218103 (2008).
  * Tailleur  _et al._ [2008] J. Tailleur, J. Kurchan, and V. Lecomte, Mapping out-of-equilibrium into equilibrium in one-dimensional transport models, Journal of Physics A: Mathematical and Theoretical 41, 505001 (2008).
  * Bertini  _et al._ [2015] L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, Macroscopic fluctuation theory, Reviews of Modern Physics 87, 593 (2015).
  * Baek  _et al._ [2018] Y. Baek, Y. Kafri, and V. Lecomte, Dynamical phase transitions in the current distribution of driven diffusive channels, Journal of Physics A: Mathematical and Theoretical 51, 105001 (2018).
  * Zakine and Vanden-Eijnden [2022] R. Zakine and E. Vanden-Eijnden, Minimum action method for nonequilibrium phase transitions, arXiv preprint arXiv:2202.06936 (2022).
  * Wu and Wang [2013] W. Wu and J. Wang, Potential and flux field landscape theory. i. global stability and dynamics of spatially dependent non-equilibrium systems, The Journal of chemical physics 139 (2013).
  * Wu and Wang [2014] W. Wu and J. Wang, Potential and flux field landscape theory. ii. non-equilibrium thermodynamics of spatially inhomogeneous stochastic dynamical systems, The Journal of Chemical Physics 141 (2014).
  * Solon  _et al._ [2018a] A. P. Solon, J. Stenhammar, M. E. Cates, Y. Kafri, and J. Tailleur, Generalized thermodynamics of phase equilibria in scalar active matter, Phys. Rev. E 97, 020602 (2018a).
  * Cates [2019] M. E. Cates, Active Field Theories, arXiv:1904.01330 [cond-mat] (2019), arXiv: 1904.01330.
  * Solon  _et al._ [2018b] A. P. Solon, J. Stenhammar, M. E. Cates, Y. Kafri, and J. Tailleur, Generalized thermodynamics of motility-induced phase separation: phase equilibria, Laplace pressure, and change of ensembles, New Journal of Physics 20, 075001 (2018b).
  * Cates  _et al._ [2010] M. E. Cates, D. Marenduzzo, I. Pagonabarraga, and J. Tailleur, Arrested phase separation in reproducing bacteria creates a generic route to pattern formation, Proceedings of the National Academy of Sciences 107, 11715 (2010).
  * Cates and Tailleur [2015] M. E. Cates and J. Tailleur, Motility-Induced Phase Separation, Annual Review of Condensed Matter Physics 6, 219 (2015).
  * Helfrich [1973] W. Helfrich, Elastic Properties of Lipid Bilayers: Theory and Possible Experiments, Zeitschrift für Naturforschung C 28, 693 (1973).
  * Bitbol  _et al._ [2012] A.-F. Bitbol, D. Constantin, and J.-B. Fournier, Bilayer Elasticity at the Nanoscale: The Need for New Terms, PLoS ONE 7, e48306 (2012).
  * Seara  _et al._ [2023] D. S. Seara, J. Colen, M. Fruchart, Y. Avni, D. Martin, and V. Vitelli, Sociohydrodynamics: data-driven modelling of social behavior, arXiv preprint arXiv:2312.17627 (2023).
  * Farley  _et al._ [1993] R. Farley, C. Steeh, T. Jackson, M. Krysan, and K. Reeves, Continued Racial Residential Segregation in Detroit: "Chocolate City, Vanilla Suburbs" Revisited, Journal of Housing Research 4, 1 (1993), publisher: American Real Estate Society.
  * Clark [2002] W. A. V. Clark, Ethnic Preferences and Ethnic Perceptions in Multi-Ethnic Settings, Urban Geography 23, 237 (2002).
  * Bobo and Zubrinsky [1996] L. Bobo and C. L. Zubrinsky, Attitudes on Residential Integration: Perceived Status Differences, Mere In-Group Preference, or Racial Prejudice?*, Social Forces 74, 883 (1996).
  * Saha  _et al._ [2020] S. Saha, J. Agudo-Canalejo, and R. Golestanian, Scalar Active Mixtures: The Nonreciprocal Cahn-Hilliard Model, Physical Review X 10, 041009 (2020).
  * Dinelli  _et al._ [2022] A. Dinelli, J. O’Byrne, A. Curatolo, Y. Zhao, P. Sollich, and J. Tailleur, Non-reciprocity across scales in active mixtures, arXiv preprint arXiv:2203.07757 (2022).
  * [76] https://github.com/rzakine/sakoda-schelling.
  * Gabaix and Koijen [2021] X. Gabaix and R. S. Koijen,  _In search of the origins of financial fluctuations: The inelastic markets hypothesis_ , Tech. Rep. (National Bureau of Economic Research, 2021).
  * Gauvin  _et al._ [2013] L. Gauvin, A. Vignes, and J.-P. Nadal, Modeling urban housing market dynamics: Can the socio-spatial segregation preserve some social diversity?, Journal of Economic Dynamics and Control 37, 1300 (2013).
  * Becharat  _et al._ [2023] A.-C. Becharat _et al._ , In preparation (2023).
  * Blume [1993] L. E. Blume, The statistical mechanics of strategic interaction, Games and Economic Behavior 5, 10.1006/game.1993.1023 (1993).
  * Brock [1993] W. A. Brock, Pathways to randomness in the economy: emergent nonlinearity and chaos in economics and finance, Estudios Economicos , 3 (1993).
  * Borghesi and Bouchaud [2007] C. Borghesi and J.-P. Bouchaud, Of songs and men: a model for multiple choice with herding, Quality & quantity 41, 557 (2007).
  * Moran  _et al._ [2020] J. Moran, A. Fosset, D. Luzzati, J.-P. Bouchaud, and M. Benzaquen, By force of habit: Self-trapping in a dynamical utility landscape, Chaos: An Interdisciplinary Journal of Nonlinear Science 30 (2020).