# Effects of odderon spin on helicity amplitudes in $pp$ elastic scattering

Prin Sawasdipol<EMAIL_ADDRESS>Khon Kaen Particle Physics and
Cosmology Theory Group (KKPaCT), Department of Physics, Faculty of Science,
Khon Kaen University,123 Mitraphap Rd., Khon Kaen, 40002, Thailand  Jingle B.
Magallanes<EMAIL_ADDRESS>Department of Physics, Mindanao
State University - Iligan Institute of Technology, Iligan City, 9200,
Philippines  Chakrit Pongkitivanichkul<EMAIL_ADDRESS>Khon Kaen Particle
Physics and Cosmology Theory Group (KKPaCT), Department of Physics, Faculty of
Science, Khon Kaen University,123 Mitraphap Rd., Khon Kaen, 40002, Thailand
Daris Samart<EMAIL_ADDRESS>corresponding author Khon Kaen Particle
Physics and Cosmology Theory Group (KKPaCT), Department of Physics, Faculty of
Science, Khon Kaen University,123 Mitraphap Rd., Khon Kaen, 40002, Thailand

###### Abstract

In recent years, the discovery of the odderon, a colorless $C$-odd gluonic
compound, has been confirmed in the TOTEM and D0 collaborations. However, the
spin quantum number of the odderon remains unidentified. In this work, we aim
to attribute a spin of $J=3$ to the odderon in $pp$ elastic scattering by
calculating the helicity amplitudes and the corresponding complex parameter
$r_{5}$, the ratio of helicity’s single-flip to non-flip amplitudes, for the
spin-3 tensor odderon with the standard spin-2 tensor pomeron exchanges. Then,
we apply these results to the constraints obtained from the STAR experiment at
RHIC. By comparing to the contributions of the spin-1 vector odderon and
spin-2 tensor pomeron, we demonstrate that the spin-3 tensor odderon, i.e.
$J=3$, provides a better explanation for the observable in $pp$ elastic
scattering.

##  I Introduction

In soft high-energy collisions, specifically hadron-hadron elastic scattering
with a large center-of-mass energy, $\sqrt{s}$, and small momentum transfer,
$\sqrt{|t|}$, the Regge theory of strong interaction phenomenology has been
widely used to explain these processes. The exchange of particles mediated by
the pomeron trajectory, as predicted by Regge theory, has proven to be
successful in describing such interactions Regge (1959). The pomeron, which
corresponds to vacuum quantum number states of the reggeon, accurately
accounts for the rising behavior observed in the total cross-section of high-
energy hadronic collisions Collins (2009); Collins and Martin (1984); Gribov
(2008). Additionally, it satisfies the Froissart bound Froissart (1961);
Lukaszuk and Martin (1967), which imposes an upper limit on the total cross-
section of two-particle scattering at high energy. The Froissart bound arises
from the Regge trajectory’s vacuum quantum numbers. Furthermore, the Regge
trajectories often yield a Chew–Frautschi plot Chew and Frautschi (1962),
which reveals the relationship between mass and spin. This significant plot is
utilized to effectively identify higher spin Regge trajectories.

Nevertheless, when considering $pp$ and $p\bar{p}$ scattering processes, a
model based solely on pomeron exchange is found to be inconsistent with the
available data. The emergence of Quantum Chromodynamics (QCD) provided
theorists with a framework to propose the inclusion of odderon exchange
alongside pomeron exchange, aiming to explain the observed rising behavior in
the total cross section for $pp$ scattering Donnachie et al. (2002); Barone
and Predazzi (2002).

The odderon was initially predicted in 1973 as a composite state of three
reggeized gluons, based on the asymptotic theorem proposed by L. Lukaszuk and
B. Nicolescu Lukaszuk and Nicolescu (1973). In detail, the odderon contributes
to scattering amplitude by decomposing the amplitude into parts associated
with charge-conjugation exchanges, specifically $C=+1$ for the pomeron and
$C=-1$ for the odderon. Two decades later, the theory of the odderon was
further developed within perturbative QCD, describing it as a bound state of
three gluons Bartels et al. (2000). Subsequently, the pomeron and odderon have
been studied through various approaches, including phenomenological Regge
theory Ewerz et al. (2014, 2016); Covolan et al. (1996); Block and Halzen
(2012); Szanyi et al. (2020); Jenkovszky (2020); Broniowski et al. (2018);
Csorgo and Szanyi (2021); Csörgő et al. (2021, 2019); Ster et al. (2015);
Block and Cahn (1985); Khoze et al. (2018), QCD-inspired models Halzen et al.
(1993); Donnachie and Landshoff (1983); Ma et al. (2001); Hu et al. (2002); He
et al. (2003); Hu et al. (2008); Zhou et al. (2006); Lu et al. (2020); Bartels
et al. (2000), and string theory or AdS/CFT correspondence Domokos et al.
(2009, 2010); Avsar et al. (2010); Hu et al. (2018); Xie et al. (2019);
Burikham and Samart (2019); Liu et al. (2022a, b); Ballon-Bayona et al.
(2017); Iatrakis et al. (2016).

In 2021, the TOTEM and D0 collaborations made a significant discovery
regarding the existence of odderon exchange in $t$-channel $pp$ and $p\bar{p}$
scattering at high energies, with a combined significance exceeding 5$\sigma$
Abazov et al. (2021). The analysis, based on a model-independent approach,
compared the $p\bar{p}$ elastic cross section at $\sqrt{s}=$ 1.96 TeV, as
measured by the D0 collaboration, with the extrapolated values obtained from
$pp$ scattering at $\sqrt{s}=$ 2.76, 7, 8, and 13 TeV, as measured by the
TOTEM collaboration. This comparison confirmed the presence of odderon
exchange. Subsequently, a theoretical study, employing model-dependent
methods, further increased the statistical significance to at least
7.08$\sigma$ by comparing the differential cross sections of elastic $pp$
scattering with the extrapolated $p\bar{p}$ scattering at LHC energies,
specifically $\sqrt{s}=$ 2.76 and 7 TeV Csorgo and Szanyi (2021).

Furthermore, it is crucial to investigate the spin structure of the pomeron
and odderon. This is because, in elastic scattering, the evaluation of the
scattering cross section requires determining the spin-dependent helicity
amplitudes of the model, using unitarity and optical theorem. Indeed, a
sophisticated framework for understanding the spin structure involves
examining the complex parameter $r_{5}$, namely the ratio of spin-flip to non-
flip amplitudes Goldberger et al. (1960); Buttimore et al. (1999). A decade
ago, the STAR experiment at RHIC conducted a study on single-spin asymmetry
measurements and derived a strong constraint on $r_{5}$ as the central value
of $(0.0017+i\,0.007)$ with its statistical and statistical+systematic
uncertainties at $\sqrt{s}=200$ GeV Adamczyk et al. (2013). Therefore, to
assess the viability of pomeron and odderon models, any spin-dependent
observable from these models that are compatible with the constraint on
$r_{5}$ can be regarded as a candidate model for pomeron and odderon
contributions.

Studies on the spin structure of the odderon date back two decades when the
relevant spin asymmetry $A_{NN}$ were demonstrated Leader and Trueman (2000).
In Ewerz et al. (2016), the possibility of pomeron spin contributions on the
$r_{5}$ parameter has been studied based on helicity amplitudes with the
Donnachie-Landshoff (DL) pomeron ansatz. Using experimental data from the STAR
at RHIC Adamczyk et al. (2013), they show that the spin-2 tensor pomeron gives
the most consistent value. Given the recent discovery of the odderon, the
possibility of odderon contribution to the $r_{5}$ parameter has been
tantalizing and become an active field of research. For example, the study on
odderon contribution based on diquark model and QCD-inspired model can be
found in Szymanowski and Zhou (2016) and Hagiwara et al. (2020) respectively.
Our recent approach on odderon as a Regge oddball spin-3 with the DL
propagator gives a theoretical basis for the odderon calculation and satisfies
several experimental results in $pp$ and $p\bar{p}$ elastic scattering
Magallanes et al. (2022). Therefore in this paper, we provide an attempt to
attribute the spin of the odderon as $J=3$ on helicity amplitudes in $pp$
elastic scattering by studying the combined exchange of pomeron and odderon
based on the aspect of spin-dependence reaction.

The objectives of this study are as follows: (1) to calculate the helicity
amplitudes of the pomeron and odderon based on Effective Field Theory (EFT)
using DL ansatz for the pomeron and odderon propagator forms, (2) to calculate
the corresponding complex parameter $r_{5}$ and then apply to the constraints
of the STAR experiment at RHIC, and (3) to evaluate the possibility of
attributing $J=3$ to the odderon that contributed with pomeron to the $pp$
elastic scattering.

We organize this paper as follows: Section II provides comprehensive details
of the theoretical frameworks for the _combined model_ of spin-2 tensor
pomeron and spin-3 tensor odderon. This model serves as the main study of our
research. In addition, the _combined model_ of spin-2 tensor pomeron and
spin-1 tensor odderon is also described and is included to complement the
evaluation of our aim. In Section III, we present the results and engage in
detailed discussions regarding the helicity amplitudes and the evaluation of
the observable $r_{5}$ for both combined models. We also compare our findings
to previous studies and assess the implications of our results. We finally
summarize our findings and draw conclusions in the Section IV highlighting the
contributions and potential future directions of research in this field.

##  II Theoretical Framework

The transition amplitude of a scattering process is of great importance as it
generally relates Quantum Field Theory (QFT) to observable in experiment.
However, obtaining unpolarized cross section for any process requires squaring
the amplitude, which leads to an exponential increase in the number of terms
and renders the calculation of complex processes infeasible.

In contrast, the helicity approach directly calculates the transition
amplitude by fixing the basis for the polarization states of external
particles. This framework incorporates a dispersion approach for crossing
relations Goldberger et al. (1957), a partial wave decomposition of the
scattering amplitudes Jacob and Wick (1959), and a connection of the
dispersion relation to partial wave amplitudes Mandelstam (1958). Together,
these components form a complete theory of nucleon-nucleon elastic scattering
for low energy, which is the momentum transfer in this work, denoted as
$\sqrt{|t|}$.

In this work, we considered the $pp$ elastic scattering process where $1,2$
are incoming and $3,4$ are outgoing particles. Each particle has its
corresponding masses $m_{i}$ four-momenta $p_{i}$ and helicities $h_{i}$,
where $i=1,..,4$ are the number of particles. According to our motivation in
this work, moreover, the contribution of the pomeron and odderon exchanges are
considered. The Feynman diagrams of the process are shown in the Fig. 1.

Figure 1:  Feynman diagram of $t$-channel $1+2\rightarrow 3+4$ elastic $pp$
scattering with the intermediate exchanges of pomeron $\mathbb{P}$ and odderon
$\mathbb{O}$.

The kinematics of the process are typically defined as

$\displaystyle\begin{aligned}
s&=(p_{1}+p_{2})^{2}=(p_{3}+p_{4})^{2}=4(\mathbf{p}^{2}+m^{2}),\\\
t&=(p_{1}-p_{3})^{2}=(p_{2}-p_{4})^{2}=-2\mathbf{p}^{2}(1-\cos(\theta)),\\\
u&=(p_{1}-p_{4})^{2}=(p_{2}-p_{3})^{2}=-2\mathbf{p}^{2}(1+\cos(\theta)),\end{aligned}$
(1)

where $s+t+u=\sum_{i=1}^{4}m_{i}^{2}$ and $p=(E,\mathbf{p})$, where $E$ is its
energy and $p$ its three-momentum. Additionally, for equal mass $m_{i}=m$ and
c.m. frame $\mathbf{p}_{1}=\mathbf{p}_{3}=\mathbf{p}$ where
$-1\leq\cos\theta\leq 1$ is the angle between the three-momenta of particles 1
and 3 in the c.m. frame, the physical region for $t$-channel scattering is
then given by $t\geq 4m^{2},~{}u\leq 0,~{}s\leq 0$ Donnachie et al. (2002).

First, we constructed the _combined spin-3 tensor-odderon_ model by employing
the contributions of the spin-2 tensor pomeron from Ewerz et al. (2016) and
spin-3 tensor odderon from Magallanes et al. (2022). The Lagrangian given from
those works read

$\displaystyle\mathcal{L}_{\mathbb{P}_{T}}(x)$
$\displaystyle=-i\,g_{\mathbb{P}_{T}}\bar{\psi}_{p}(x)\,\mathcal{G}_{(\mu\nu)}^{\mu^{\prime}\nu^{\prime}}\,\gamma_{\mu^{\prime}}\,\overleftrightarrow{\partial_{\nu^{\prime}}}\,\psi_{p}(x)\,\mathbb{P}^{\mu\nu}(x),$
(2) $\displaystyle\mathcal{L}_{\mathbb{O}_{T}}(x)$
$\displaystyle=-\frac{g_{\mathbb{O}_{T}}}{{M}_{0}}\,\mathcal{G}_{(\mu\nu\rho)}^{\mu^{\prime}\nu^{\prime}\rho^{\prime}}\,\bar{\psi}_{p}(x)\,\gamma_{\mu^{\prime}}\,\overleftrightarrow{\partial_{\nu^{\prime}}}\,\overleftrightarrow{\partial_{\rho^{\prime}}}\,\psi_{p}(x)\,\mathbb{O}^{\mu\nu\rho}(x),$
(3)

where $g_{\mathbb{P}_{T}}$ and $g_{\mathbb{O}_{T}}$ are coupling constants and
they carry the mass dimension as [mass]-1, $\mathbb{P}^{\mu\nu}(x)$ and
$\mathbb{O}^{\mu\nu\rho}(x)$ denote its descriptions, and $\mathbb{P}_{T}$ and
$\mathbb{O}_{T}$ denote the spin-2 tensor pomeron and spin-3 tensor odderon,
respectively. The mass parameter $M_{0}$ is a free parameter and is introduced
for correcting mass dimension. $\psi_{p}$ is the Dirac field operator of the
proton.
$\overleftrightarrow{\partial_{\mu}}\equiv\overleftarrow{\partial_{\mu}}-\overrightarrow{\partial_{\mu}}$.
The totally symmetric operators given by

$\displaystyle\mathcal{G}_{(\mu\nu)}^{\mu^{\prime}\nu^{\prime}}$
$\displaystyle=\frac{1}{2!}\left(g_{\mu}^{\mu^{\prime}}g_{\nu}^{\nu^{\prime}}+g_{\nu}^{\mu^{\prime}}g_{\mu}^{\nu^{\prime}}\right),$
(4)
$\displaystyle\mathcal{G}_{(\mu\nu\rho)}^{\mu^{\prime}\nu^{\prime}\rho^{\prime}}$
$\displaystyle=\frac{1}{3!}\left(g_{\mu}^{\mu^{\prime}}g_{\nu}^{\nu^{\prime}}g_{\rho}^{\rho^{\prime}}+g_{\nu}^{\mu^{\prime}}g_{\rho}^{\nu^{\prime}}g_{\mu}^{\rho^{\prime}}+g_{\rho}^{\mu^{\prime}}g_{\mu}^{\nu^{\prime}}g_{\nu}^{\rho^{\prime}}+g_{\nu}^{\mu^{\prime}}g_{\mu}^{\nu^{\prime}}g_{\rho}^{\rho^{\prime}}+g_{\mu}^{\mu^{\prime}}g_{\rho}^{\nu^{\prime}}g_{\nu}^{\rho^{\prime}}+g_{\rho}^{\mu^{\prime}}g_{\nu}^{\nu^{\prime}}g_{\mu}^{\rho^{\prime}}\right),$
(5)

are used to ensure the totally-symmetric Lorentz indices, $(\mu\nu)$ and
$(\mu\nu\rho)$, respectively. Note that the metric tensor signature in this
work is $g_{\mu\nu}={\rm diag}(1,-1,-1,-1)$.

With using of the standard QFT method Peskin and Schroeder (1995), the
vertices of the interaction Lagrangian in Eqs.(2) and (3) are formulated as

$\displaystyle i\Gamma_{\mu\nu}^{\mathbb{P}_{T}}(p^{\prime},p)$
$\displaystyle=-i\,g_{\mathbb{P}_{T}}\,\mathcal{G}_{(\mu\nu)}^{\mu^{\prime}\nu^{\prime}}\,\gamma_{\mu^{\prime}}\,(p^{\prime}+p)_{\nu^{\prime}}\,\mathcal{F}_{\mathbb{P}_{T}}\left[(p^{\prime}-p)^{2}\right],$
(6) $\displaystyle i\Gamma_{\mu\nu\rho}^{\mathbb{O}_{T}}(p^{\prime},p)$
$\displaystyle=-i\,\frac{g_{\mathbb{O}_{T}}}{{M}_{0}}\,\mathcal{G}_{(\mu\nu\rho)}^{\mu^{\prime}\nu^{\prime}\rho^{\prime}}\,\gamma_{\mu^{\prime}}\,(p^{\prime}+p)_{\nu^{\prime}}\,(p^{\prime}+p)_{\rho^{\prime}}\,\mathcal{F}_{\mathbb{O}_{T}}\left[(p^{\prime}-p)^{2}\right],$
(7)

where $p$ and $p^{\prime}$ are the momentum of incoming and outgoing particles
respectively. Also, the coupling constants
$g_{\mathbb{P}_{T}}=3\beta_{\mathbb{P}NN}=3\times 1.87$ GeV-1 and
$g_{\mathbb{O}_{T}}=13.80$ GeV-1 are obtained from Ewerz et al. (2014) and
Magallanes et al. (2022), respectively. The pomeron and odderon form factors
are given by

$\displaystyle\mathcal{F}_{\mathbb{P}_{T}}(t)$
$\displaystyle=\left(1-\frac{t}{4\,m_{p}^{2}}\frac{\mu_{p}}{\mu_{N}}\right)\left(1-\frac{t}{4\,m_{p}^{2}}\right)^{-1}\left(1-\frac{t}{m_{D}^{2}}\right)^{-2},$
(8) $\displaystyle\mathcal{F}_{\mathbb{O}_{T}}(t)$
$\displaystyle=\left(1-\frac{A\,t}{4\,m_{p}^{2}}\frac{\mu_{p}}{\mu_{N}}\right)\left(1-\frac{B\,t}{4\,m_{p}^{2}}\right)^{-1}\left(1-\frac{C\,t}{m_{D}^{2}}\right)^{-2},$
(9)

with their corresponding normalization of $\mathcal{F}_{\mathbb{P}_{T}}(0)=1$
and $\mathcal{F}_{\mathbb{O}_{T}}(0)=1$, where the proton mass $m_{p}=0.978$
GeV, the ratio of magnetic moment of proton to nucleon
$\mu_{p}/\mu_{N}=2.7928$, and the squared dipole mass $m_{D}^{2}=0.71$ GeV2
are given from Ewerz et al. (2014). Note that the form factor of the pomeron
has no free parameters whereas for odderon the form factor is modified by
three parameters: $A=-1.671$, $B=-1.855$, and $C=0.493$ which is obtained from
fitting the model to several experiments, as shown in Magallanes et al.
(2022).

The corresponding propagators read

$\displaystyle i\Delta_{\mu\nu,\rho\sigma}^{\mathbb{P}_{T}}(s,t)$
$\displaystyle=\frac{1}{4\,s}\left[g_{\mu\rho}\,g_{\nu\sigma}+g_{\mu\sigma}\,g_{\nu\rho}-\frac{1}{2}g_{\mu\nu}\,g_{\rho\sigma}\right](-i\,s\,\alpha^{\prime}_{\mathbb{P}_{T}})^{\alpha_{\mathbb{P}_{T}}(t)-1},$
(10) $\displaystyle
i\Delta_{\mu\nu\lambda,\rho\sigma\tau}^{\mathbb{O}_{T}}(s,t)$
$\displaystyle=-i\,\frac{{M}_{0}^{2}}{6\,s^{2}}\left[\sum_{C}g_{\mu\rho}(g_{\mu\sigma}\,g_{\lambda\tau}+g_{\nu\tau}\,g_{\lambda\sigma})-\frac{1}{2}\sum_{C}g_{\mu\nu}\,g_{\lambda\rho}\,g_{\sigma\tau}\right](-i\,s\,\alpha^{\prime}_{\mathbb{O}_{T}})^{\alpha_{\mathbb{O}_{T}}(t)-1},$
(11)

where $\epsilon$ and $\alpha^{\prime}$ denote the interception deviation and
trajectory slope, respectively. For spin-2 tensor pomeron, the trajectory is
given by
$\alpha_{\mathbb{P}_{T}}(t)=1+\epsilon_{\mathbb{P}_{T}}+\alpha_{\mathbb{P}_{T}}^{\prime}t$
with the $\epsilon_{\mathbb{P}_{T}}=0.0808$ and
$\alpha_{\mathbb{P}_{T}}^{\prime}=0.25$ GeV2 Ewerz et al. (2014). For spin-3
tensor odderon, the trajectory is given by
$\alpha_{\mathbb{O}_{T}}(t)=1+\epsilon_{\mathbb{O}_{T}}+\alpha_{\mathbb{O}_{T}}^{\prime}t$
with the $\epsilon_{\mathbb{O}_{T}}=0.062$ and
$\alpha_{\mathbb{O}_{T}}^{\prime}=0.189$ GeV2 Magallanes et al. (2022).

Next, we will calculate the corresponding helicity amplitudes in the following
forms,

$\displaystyle\begin{aligned} i\langle
h_{3}h_{4}|\mathcal{T}_{\mathbb{P}_{T}}|h_{1}h_{2}\rangle&=\langle
p(p_{3},h_{3})p(p_{4},h_{4})|:\mathbb{T}\exp\left(i\int\mathcal{L}_{\mathbb{P}_{T}}(x)d^{4}x\right):|p(p_{1},h_{1})p(p_{2},h_{2})\rangle\\\
&=\bar{u}(p_{3},h_{3})\,i\Gamma_{\mu\nu}^{\mathbb{P}_{T}}(p_{3},p_{1})\,u(p_{1},h_{1})\,i\Delta_{\mathbb{P}_{T}}^{\mu\nu,\rho\sigma}(s,t)\,\bar{u}(p_{4},h_{4})\,i\Gamma_{\rho\sigma}^{\mathbb{P}_{T}}(p_{4},p_{2})\,u(p_{2},h_{2}),\end{aligned}$
(12) $\displaystyle\begin{aligned} i\langle
h_{3}h_{4}|\mathcal{T}_{\mathbb{O}_{T}}|h_{1}h_{2}\rangle&=\langle
p(p_{3},h_{3})p(p_{4},h_{4})|:\mathbb{T}\exp\left(i\int\mathcal{L}_{\mathbb{O}_{T}}(x)d^{4}x\right):|p(p_{1},h_{1})p(p_{2},h_{2})\rangle\\\
&=\bar{u}(p_{3},h_{3})\,i\Gamma_{\mu\nu\lambda}^{\mathbb{O}_{T}}(p_{3},p_{1})\,u(p_{1},h_{1})\,i\Delta_{\mathbb{O}_{T}}^{\mu\nu\lambda,\rho\sigma\tau}(s,t)\,\bar{u}(p_{4},h_{4})\,i\Gamma_{\rho\sigma\tau}^{\mathbb{O}_{T}}(p_{4},p_{2})\,u(p_{2},h_{2}),\end{aligned}$
(13)

where $u(p,h)$ is the Dirac helicity spinor defined in the Appendix A,
$\mathbb{T}$ is the time ordered operator, and $:\;:$ is the normal order
operator in the standard QFT Peskin and Schroeder (1995).

According to the invariant under the parity conservation ($P$), charge
conjugation ($C$), time reversal ($T$), and conservation of total spin, total
helicity amplitudes of the two-body fermion-fermion scattering are reduced
from sixteen down to five independent amplitudes, as seen in Donnachie et al.
(2002). Thus, the scattering amplitudes of the odderon are included, with a
minus sign, in the exchange as shown in Magallanes et al. (2022). The
scattering amplitudes of the $pp$ and $p\bar{p}$ are written by

$\displaystyle\langle h_{3}h_{4}|\mathcal{T}_{pp}|h_{1}h_{2}\rangle$
$\displaystyle=$ $\displaystyle\langle
h_{3}h_{4}|\mathcal{T}_{\mathbb{P}}-\mathcal{T}_{\mathbb{O}}|h_{1}h_{2}\rangle\,,$
(14) $\displaystyle\langle
h_{3}h_{4}|\mathcal{T}_{p\bar{p}}|h_{1}h_{2}\rangle$ $\displaystyle=$
$\displaystyle\langle
h_{3}h_{4}|\mathcal{T}_{\mathbb{P}}+\mathcal{T}_{\mathbb{O}}|h_{1}h_{2}\rangle\,,$
(15)

where $\mathcal{T}_{\mathbb{P}}$ and $\mathcal{T}_{\mathbb{O}}$ are the
transition matrices of the pomeron and odderon, respectively. Consequently,
the generic form of the independent helicity amplitudes defined in Goldberger
et al. (1960); Buttimore et al. (1999); Donnachie et al. (2002) as
$\phi_{i}=\langle h_{3}h_{4}|\mathcal{T}|h_{1}h_{2}\rangle$ are then
decomposed into pomeron and odderon contributions related to the transition
matrix $\mathcal{T}$ as given by,

$\displaystyle\begin{aligned}
\phi_{1}^{TT}(s,t)=\langle++|\mathcal{T}_{pp}^{TT}|++\rangle=\langle++|\mathcal{T}_{\mathbb{P}_{T}}-\mathcal{T}_{\mathbb{O}_{T}}|++\rangle,\\\
\phi_{2}^{TT}(s,t)=\langle++|\mathcal{T}_{pp}^{TT}|--\rangle=\langle++|\mathcal{T}_{\mathbb{P}_{T}}-\mathcal{T}_{\mathbb{O}_{T}}|--\rangle,\\\
\phi_{3}^{TT}(s,t)=\langle+-|\mathcal{T}_{pp}^{TT}|+-\rangle=\langle+-|\mathcal{T}_{\mathbb{P}_{T}}-\mathcal{T}_{\mathbb{O}_{T}}|+-\rangle,\\\
\phi_{4}^{TT}(s,t)=\langle+-|\mathcal{T}_{pp}^{TT}|-+\rangle=\langle+-|\mathcal{T}_{\mathbb{P}_{T}}-\mathcal{T}_{\mathbb{O}_{T}}|-+\rangle,\\\
\phi_{5}^{TT}(s,t)=\langle++|\mathcal{T}_{pp}^{TT}|+-\rangle=\langle++|\mathcal{T}_{\mathbb{P}_{T}}-\mathcal{T}_{\mathbb{O}_{T}}|+-\rangle,\end{aligned}$
(16)

where plus ($+$) and minus ($-$) signs in the states correspond to
right-$(\chi_{\uparrow})$ and left-$(\chi_{\downarrow})$ handed helicities,
mentioned in the Appendix A, respectively. Physically, $\phi_{1}$ and
$\phi_{3}$ are spin-flip amplitudes, $\phi_{2}$ and $\phi_{4}$ are double-flip
amplitudes, and $\phi_{5}$ is single-flip amplitude. The superscript “$TT$”
denotes the consecutive letters as tensors for spin-2 tensor pomeron and
spin-3 tensor odderon. Note that we will apply the Regge limit $(s\gg t,m)$ to
each amplitude.

Finally, the general form of the ratio of helicity’s spin-flip to non-flip
amplitudes given by Buttimore et al. (1999), see also Ewerz et al. (2016), as

$\displaystyle r_{5}^{TT}(s,t)$
$\displaystyle=\frac{2\,m_{p}\,\phi_{5}^{TT}(s,t)}{\sqrt{-t}\,\textrm{Im}[\phi_{1}^{TT}(s,t)+\phi_{3}^{TT}(s,t)]}.$
(17)

In addition, to get a convenient form for presenting the results of this work
and comparing them to others, we defined the reduced amplitudes as

$\displaystyle\begin{aligned}
\hat{\phi}_{i}^{\mathbb{P}_{T}}(s,t)&=\phi_{i}^{\mathbb{P}_{T}}(s,t)/\mathcal{C}_{\mathbb{P}_{T}}(s,t),\\\
\hat{\phi}_{i}^{\mathbb{O}_{T}}(s,t)&=\phi_{i}^{\mathbb{O}_{T}}(s,t)/\mathcal{C}_{\mathbb{O}_{T}}(s,t),\end{aligned}$
(18)

where $i=1,..,5$ is the index of independent helicity amplitudes as used in
(16), and the common factors are defined as

$\displaystyle\begin{aligned}
\mathcal{C}_{\mathbb{P}_{T}}(s,t)&=i\,g_{\mathbb{P}_{T}}^{2}\mathcal{F}_{\mathbb{P}_{T}}^{2}(t)\frac{1}{4\,s}(-i\,s\,\alpha^{\prime}_{\mathbb{P}_{T}})^{\alpha_{\mathbb{P}_{T}}(t)-1},\\\
\mathcal{C}_{\mathbb{O}_{T}}(s,t)&=\frac{8}{9}\,g_{\mathbb{O}_{T}}^{2}\mathcal{F}_{\mathbb{O}_{T}}^{2}(t)\frac{1}{6\,s}(-i\,s\,\alpha^{\prime}_{\mathbb{O}_{T}})^{\alpha_{\mathbb{O}_{T}}(t)-1}.\end{aligned}$
(19)

For the odderon as a spin-1 vector particle, we, therefore, constructed the
_combined spin-1 vector odderon_ model in order to calculate and compare the
results to that of the _combined spin-3 tensor odderon_.

The Lagrangian of the spin-1 vector odderon from Ewerz et al. (2014) reads

$\displaystyle\mathcal{L}_{\mathbb{O}_{V}}(x)$
$\displaystyle=-g_{\mathbb{O}_{V}}M_{0}\,\bar{\psi}_{p}(x)\,\gamma_{\mu}\psi_{p}(x)\,\mathbb{O}^{\mu}(x),$
(20)

where $g_{\mathbb{O}_{V}}=$ is the coupling constant, $\mathbb{O}^{\mu}(x)$ is
the description of the spin-1 vector odderon. The corresponding vertex reads

$\displaystyle i\Gamma_{\mu}^{\mathbb{O}_{V}}(p^{\prime},p)$
$\displaystyle=-i\,g_{\mathbb{O}_{V}}\,M_{0}\,\gamma_{\mu}\,\mathcal{F}_{\mathbb{O}_{V}}\left[(p^{\prime}-p)^{2}\right],$
(21)

where the form factor is defined the same as the pomeron,

$\displaystyle\mathcal{F}_{\mathbb{O}_{V}}(t)$
$\displaystyle=\left(1-\frac{t}{4m_{p}^{2}}\frac{\mu_{p}}{\mu_{N}}\right)\left(1-\frac{t}{4m_{p}^{2}}\right)^{-1}\left(1-\frac{t}{m_{D}^{2}}\right)^{-2},$
(22)

with normalization $\mathcal{F}_{\mathbb{O}_{V}}(0)=1$. The corresponding
propagator reads

$\displaystyle i\Delta_{\mu\nu}^{\mathbb{O}_{V}}(s,t)$
$\displaystyle=-i\frac{g_{\mu\nu}}{M_{0}^{2}}(-i\,s\,\alpha^{\prime}_{\mathbb{O}_{V}})^{\alpha_{\mathbb{O}_{V}}(t)-1},$
(23)

where the odderon trajectory is
$\alpha_{\mathbb{O}_{V}}(t)=1+\epsilon_{\mathbb{O}_{V}}+\alpha^{\prime}_{\mathbb{O}_{V}}t$
with both $\epsilon_{\mathbb{O}_{V}}=0.0808$ and
$\alpha^{\prime}_{\mathbb{O}_{V}}=0.25$ GeV2 are also taken from Ewerz et al.
(2014). Thus, a form of the corresponding helicity amplitudes are

$\displaystyle\begin{aligned} i\langle
h_{3}h_{4}|\mathcal{T}_{\mathbb{O}_{V}}|h_{1}h_{2}\rangle&=\langle
p(p_{3},h_{3})p(p_{4},h_{4})|:\mathbb{T}\exp\left(i\int\mathcal{L}_{\mathbb{O}_{V}}(x)d^{4}x\right):|p(p_{1},h_{1})p(p_{2},h_{2})\rangle\\\
&=\bar{u}(p_{3},h_{3})\,i\Gamma_{\mu}^{\mathbb{O}_{V}}(p_{3},p_{1})\,u(p_{1},h_{1})\,i\Delta_{\mathbb{O}_{V}}^{\mu\nu}(s,t)\,\bar{u}(p_{4},h_{4})\,i\Gamma_{\nu}^{\mathbb{O}_{V}}(p_{4},p_{2})\,u(p_{2},h_{2}).\end{aligned}$
(24)

Consequently, we defined the independent helicity amplitudes of the spin-1
vector odderon in the same manner as the previous case,

$\displaystyle\begin{aligned}
\phi_{1}^{TV}(s,t)=\langle++|\mathcal{T}_{pp}^{TV}|++\rangle=\langle++|\mathcal{T}_{\mathbb{P}_{T}}-\mathcal{T}_{\mathbb{O}_{V}}|++\rangle,\\\
\phi_{2}^{TV}(s,t)=\langle++|\mathcal{T}_{pp}^{TV}|--\rangle=\langle++|\mathcal{T}_{\mathbb{P}_{T}}-\mathcal{T}_{\mathbb{O}_{V}}|--\rangle,\\\
\phi_{3}^{TV}(s,t)=\langle+-|\mathcal{T}_{pp}^{TV}|+-\rangle=\langle+-|\mathcal{T}_{\mathbb{P}_{T}}-\mathcal{T}_{\mathbb{O}_{V}}|+-\rangle,\\\
\phi_{4}^{TV}(s,t)=\langle+-|\mathcal{T}_{pp}^{TV}|-+\rangle=\langle+-|\mathcal{T}_{\mathbb{P}_{T}}-\mathcal{T}_{\mathbb{O}_{V}}|-+\rangle,\\\
\phi_{5}^{TV}(s,t)=\langle++|\mathcal{T}_{pp}^{TV}|+-\rangle=\langle++|\mathcal{T}_{\mathbb{P}_{T}}-\mathcal{T}_{\mathbb{O}_{V}}|+-\rangle,\end{aligned}$
(25)

where the superscript “TV” represents the spin-2 tensor pomeron and the spin-1
vector odderon exchanges in $pp$ elastic scattering and the corresponding
result of the $r_{5}$ parameter reads

$\displaystyle r_{5}^{TV}(s,t)$
$\displaystyle=\frac{2\,m_{p}\,\phi_{5}^{TV}(s,t)}{\sqrt{-t}\,\textrm{Im}[\phi_{1}^{TV}(s,t)+\phi_{3}^{TV}(s,t)]}.$
(26)

In the same manner, we defined the helicity amplitudes in the following form,

$\displaystyle\hat{\phi}_{i}^{\mathbb{O}_{V}}(s,t)$
$\displaystyle=\phi_{i}^{\mathbb{O}_{V}}(s,t)/\mathcal{C}_{\mathbb{O}_{V}}(s,t),$
(27)

where the corresponding common factor reads

$\displaystyle\mathcal{C}_{\mathbb{O}_{V}}(s,t)$
$\displaystyle=g_{\mathbb{O}_{V}}^{2}\,\mathcal{F}_{\mathbb{O}_{V}}^{2}(t)\frac{1}{4\,s}(-i\,s\,\alpha^{\prime}_{\mathbb{O}_{V}})^{\alpha_{\mathbb{O}_{V}}(t)-1}.$
(28)

In this section, we have formulated the relevant observable from the helicity
amplitudes in $pp$ elastic scattering. We have also categorized the complex
parameter $r_{5}$ as the results of two models: (1) the $r_{5}^{TT}$ from the
spin-2 tensor pomeron with the spin-3 tensor odderon exchanges and (2) the
$r_{5}^{TV}$ from the spin-2 tensor pomeron with the spin-1 vector odderon
exchanges. We will compare the theoretical results from the $r_{5}^{TT}$ and
$r_{5}^{TV}$ with the experimental data in the next section.

##  III  Results and Discussion

The purpose of this research is to explore the possibility that the odderon
with the spin $J=3$, alongside the pomeron, contributes to $pp$ elastic
scattering. Therefore, the observable as a complex variable $r_{5}$, namely
the ratio of spin-flip to non-flip amplitudes, of the combined exchange of the
odderon and pomeron is evaluated from their corresponding helicity amplitudes.

First, we have calculated the helicity amplitudes of the spin-2 tensor pomeron
($\mathbb{P}_{T}$) and spin-3 tensor odderon ($\mathbb{O}_{T}$) contributions,
by using (16), in the Regge limit $(s\gg t,m)$. The results are shown in the
Table 1.

|  | $\mathbb{P}_{T}$ |  | $\mathbb{O}_{T}$ |  | $\mathbb{O}_{V}$
---|---|---|---|---|---|---
$\hat{\phi}_{1}(s,t)$ |  | $8s^{2}$ |  | $9s^{2}$ |  | $8s^{2}$
$\hat{\phi}_{2}(s,t)$ |  | $10m_{p}^{2}t$ |  | $10m_{p}^{2}t$ |  | $16m_{p}^{2}t$
$\hat{\phi}_{3}(s,t)$ |  | $8s^{2}$ |  | $9s^{2}$ |  | $8s^{2}$
$\hat{\phi}_{4}(s,t)$ |  | $-10m_{p}^{2}t$ |  | $-10m_{p}^{2}t$ |  | $-16m_{p}^{2}t$
$\hat{\phi}_{5}(s,t)$ |  | $-8m_{p}s\sqrt{-t}$ |  | $-9m_{p}s\sqrt{-t}$ |  | $-8m_{p}s\sqrt{-t}$
Table 1: The reduced helicity amplitudes of the pomeron ($\mathbb{P}_{T}$) and
odderon ($\mathbb{O}_{T}$ and $\mathbb{O}_{V}$) contributions.

According to (17), the result of the spin-3 tensor odderon is written as

$\displaystyle\begin{aligned}
r_{5}^{TT}(s,t)&=\frac{m_{p}^{2}(9\,\mathcal{C}_{\mathbb{O}_{T}}(s,t)-8\,\mathcal{C}_{\mathbb{P}_{T}}(s,t))}{s\,{\rm
Im}(8\,\mathcal{C}_{\mathbb{P}_{T}}(s,t)-9\,\mathcal{C}_{\mathbb{O}_{T}}(s,t))}\\\
&=(5.3-i\,2.2)\times 10^{-5},\end{aligned}$ (29)

where the numerical result in the last line comes from the STAR experiment at
RHIC Adamczyk et al. (2013) as $\sqrt{s}=200$ GeV and no $t$-dependent
($t=0$).

Similarly we have worked out the spin-1 vector odderon ($\mathbb{O}_{V}$)
contribution. The helicity amplitudes from (25) are calculated using the
results in the Table 1. Thus, the corresponding result according to (26) is
given by

$\displaystyle\begin{aligned}
r_{5}^{TV}(s,t)&=\frac{m_{p}^{2}(\mathcal{C}_{\mathbb{O}_{V}}(s,t)-\mathcal{C}_{\mathbb{P}_{T}}(s,t))}{s\,{\rm
Im}(\mathcal{C}_{\mathbb{P}_{T}}(s,t)-\mathcal{C}_{\mathbb{O}_{V}}(s,t))},\\\
&=(1.7-i\,2.2)\times 10^{-5}.\end{aligned}$ (30)

We then evaluated our results of $r_{5}^{TT}$ by comparing to that of the
spin-2 tensor pomeron studied in Ewerz et al. (2016). The constraints on the
$r_{5}$ from the STAR experiment at RHIC Adamczyk et al. (2013) are applied to
these results for their validity. As shown in the Fig. 2, the central value of
the complex parameter $r_{5}$ is obtained from the experimental study of spin-
single asymmetry $A_{N}$ of polarized $pp$ elastic scattering at
$\sqrt{s}=200$ GeV Adamczyk et al. (2013) as $(170.0+i\,170.0)\times 10^{-5}$
with the corresponding statistical uncertainties $\sigma_{\rm
stat}=\pm(170.0+i\,300.0)\times 10^{-5}$ and statistical+systematic
uncertainties $\sigma_{\rm stat+sys}=\pm(630.0+i\,570.0)\times 10^{-5}$.

Figure 2:  (color online) The statistical (solid ellipse and cross) and
statistical+systematic uncertainties (dashed ellipse and cross) of the $r_{5}$
Adamczyk et al. (2013). The _pure spin-2 tensor pomeron_ (black circle) Ewerz
et al. (2014) is $r_{5}^{\mathbb{P}_{T}}$. The _combined spin-1 vector
odderon_ (blue square) is $r_{5}^{TV}$ and _combined spin-3 tensor odderon_
(red diamond) is $r_{5}^{TT}$.

The value from the _pure spin-2 tensor pomeron_ model is given as a benchmark,
$r_{5}^{\mathbb{P}_{T}}=-(0.28+i\,2.2)\times 10^{-5}$ Ewerz et al. (2016).
Notice that the real part of $r_{5}^{\mathbb{P}_{T}}$ is negative which is
similar to the result from the study of QCD-inspired dipole model Hagiwara et
al. (2020). The value lies outside the range of the statistical uncertainties
while it remains inside the statistical+systematic uncertainties.

The main result of this work is the value of $r_{5}^{TT}$ of the _combined
spin-3 tensor odderon_ model which has a positive real part. This comes from
the fact that the contribution of the odderon has the opposite sign against
the pomeron, as defined in (16) and (25), see also Donnachie et al. (2002);
Magallanes et al. (2022). This value is compatible with the constraints given
from the STAR experiment at RHIC Adamczyk et al. (2013). Importantly, the
clear tendency toward the central value obtained from the experiment indicates
the improvement of the model compared to the _pure spin-2 tensor pomeron_
model Ewerz et al. (2016).

In addition, the value of $r_{5}^{TV}$ of the _combined spin-1 vector odderon_
model is properly valid the same order as in the _pure spin-2 tensor pomeron_
and _combined spin-3 tensor odderon_ models. Superficially, the real part has
the positive value as the previous case. However, despite the tendency, its
distance to the central value obtained from the experiment is farther than the
result of the _combined spin-3 tensor odderon_ model. This means that the spin
of odderon is likely to be $J=3$ rather than $J=1$.

##  IV Conclusion

In spite of the discovery, the spin-dependence aspect of the odderon
contribution in $pp$ elastic scattering has not yet been thoroughly
investigated. Due to the fact that the odderon is respected as a charge-
conjugation partner of pomeron, the observable of the model naturally requires
the inclusion of both pomeron and odderon contributions.

Therefore, we have constructed the combined models of spin-3 tensor odderon
and spin-1 vector odderon with the spin-2 tensor pomeron by using the DL
ansatz for establishing the propagator forms. In each model, the complex
parameter $r_{5}$, as a result, has been calculated from the helicity
amplitudes of the models. Both models have been applied to the constraints
given from the STAR experiment at RHIC Adamczyk et al. (2013). In detail, at
first, we have reproduced the parameter $r_{5}^{\mathbb{P}_{T}}$ of the _pure
spin-2 tensor pomeron_ with the free parameters given from Ewerz et al.
(2016). The result of the spin-2 tensor pomeron is identical to that of in
Ewerz et al. (2016) and compatible with the study of Hagiwara et al. (2020).
Next, the parameter $r_{5}^{TT}$ have also been calculated from the combined
helicity amplitudes of the spin-2 tensor pomeron and the spin-3 tensor
odderon. This _combined spin-3 tensor odderon_ model is the main target of our
study. The result of the _combined spin-3 tensor odderon_ model is valid with
the constraints obtained from the experiment. Also, given that the positive
sign of the Re $r_{5}^{TT}$ against the negative sign of the Re
$r_{5}^{\mathbb{P}_{T}}$ reflects the charge conjugation partner to the
contribution of the $C$-even pomeron and $C$-odd odderon. Then, we worked out
the parameter $r_{5}^{TV}$ of the spin-1 vector odderon in the same manner as
the spin-3 tensor odderon. We have found that, despite the same tendency, the
result of the _combined spin-1 vector odderon_ model is farther from the
central value than the _combined spin-3 tensor odderon_ model. Finally, we
plotted all the results with the constraints from the experiment and drawn the
conclusion that the spin-3 tensor odderon is more likely to contribute to $pp$
elastic scattering.

The contribution of the pomeron and odderon are expected to be dominated at
TeV energy scale in $pp$ elastic scattering. In fact the parameters of our
model have been chosen from the best fit of various experiments Magallanes et
al. (2022). However, in this work the model is subjected to the constrain from
the lower energy experiment of $\sqrt{s}=200$ GeV which is the only available
data for the observable $r_{5}$ of spin-dependence $pp$ elastic scattering at
the present time. Therefore the future experiments at higher energy would
shine the light on the odderon contribution in helicity amplitudes

The fact that our results are still outside the statistical error of the
constraints given from the STAR experiment at RHIC indicates that the other
relevant exchanges could be introduced into the reaction at lower energy scale
in order to explain the experimental result. Therefore, as described in Ewerz
et al. (2014), the other contributions for the exchange in the reaction, such
as $f_{2R},a_{2R},\omega_{R},\rho_{R}$ and $\gamma$, are of interest to study
further and we leave these topics for future works.

##  Appendix A The Helicity Spinors and Kinematics

The central results in this work are helicity amplitudes of $pp$ elastic
scattering with the pomeron and odderon exchanges. Therefore, in order to make
the calculations in this work transparently and well organized, it is worth
fixing the normalization of the Dirac spinors and their helicity conventions.
With $\bar{u}_{r}(p,m)\,u_{s}(p,m)=2m\delta_{rs}$ normalization, one can write
the Dirac helicity spinors as

$\displaystyle\begin{aligned}
u_{\uparrow}(p,m)&=\sqrt{E+m}\begin{bmatrix}\chi_{\uparrow}\\\
\kappa\chi_{\uparrow}\end{bmatrix},\quad\quad
u_{\downarrow}(p,m)=\sqrt{E+m}\begin{bmatrix}\chi_{\downarrow}\\\
-\kappa\chi_{\downarrow}\end{bmatrix},\\\
v_{\uparrow}(p,m)&=\sqrt{E+m}\begin{bmatrix}-\kappa\chi_{\downarrow}\\\
\chi_{\downarrow}\end{bmatrix},\quad\quad
v_{\downarrow}(p,m)=\sqrt{E+m}\begin{bmatrix}\kappa\chi_{\uparrow}\\\
\chi_{\uparrow}\end{bmatrix},\end{aligned}$ (31)

where $\kappa=\frac{\lvert\vec{p}\rvert}{E+m}=\frac{\gamma\beta}{\gamma+1}$,
$\gamma=1/\sqrt{1-\beta^{2}}$, and $\beta=\sqrt{1-4m^{2}/s}$. While the two-
component spinors read

$\displaystyle\chi_{\uparrow}=\begin{bmatrix}\cos(\theta/2)\\\
e^{i\phi}\sin(\theta/2)\end{bmatrix},\quad\quad\chi_{\downarrow}=\begin{bmatrix}-\sin(\theta/2)\\\
e^{i\phi}\cos(\theta/2)\end{bmatrix},$ (32)

where $\phi$ is the angle in the spin-related plane of the particles. The four
momenta of the particles are read, (Figure 1)

$\displaystyle\begin{aligned} p_{1}^{\mu}&=\frac{\sqrt{s}}{2}(1,0,0,\beta),\\\
p_{2}^{\mu}&=\frac{\sqrt{s}}{2}(1,0,0,-\beta),\\\
p_{3}^{\mu}&=\frac{\sqrt{s}}{2}(1,\beta\sin{\theta},0,\beta\cos{\theta}),\\\
p_{4}^{\mu}&=\frac{\sqrt{s}}{2}(1,-\beta\sin{\theta},0,-\beta\cos{\theta}).\end{aligned}$
(33)

In addition, all momenta are fixed the directions on the c.m. frame scattering
plane as

$\displaystyle\begin{aligned} p_{1}&:\,\theta=0\,,\,\phi=0\,,\\\
p_{2}&:\,\theta=\pi\,,\,\phi=\pi\,,\\\ p_{3}&:\,\theta\,,\,\phi=0\,,\\\
p_{4}&:\,\theta=\pi-\theta\,,\,\phi=\pi.\end{aligned}$ (34)

In the Regge limit $(s\gg t,m)$, moreover, the trigonometric functions of the
c.m. frame angle can be written in terms of the Mandelstam variables as

$\displaystyle\begin{aligned} \sin({\theta/2})&\approx\sqrt{-t/s},\\\
\cos({\theta/2})&\approx 1,\\\ \sin(\theta)&\approx 2\sqrt{-{t/s}},\\\
\cos(\theta)&\approx 1,\\\ \cos(2\theta)&\approx 0.\end{aligned}$ (35)

We will use all notations and conventions introduced above carry out the
calculations throughout this work.

###### Acknowledgements.

DS is supported by the Fundamental Fund of Khon Kaen University and DS has
received funding support from the National Science, Research and Innovation
Fund. The Mindanao State University - Iligan Institute of Technology is also
acknowledged through its research and extension support extended to J.B.
Magallanes for his travel to Khon Kaen University, Thailand. PS, CP, and DS
are financially supported by the National Astronomical Research Institute of
Thailand (NARIT). CP and DS are supported by Thailand NSRF via PMU-B [grant
number B05F650021]. CP is also supported by Fundamental Fund 2565 of Khon Kaen
University and Research Grant for New Scholar, Office of the Permanent
Secretary, Ministry of Higher Education, Science, Research and Innovation
under contract no. RGNS64-043. The authors acknowledge the National Science
and Technology Development Agency, National e-Science Infrastructure
Consortium, Chulalongkorn University and the Chulalongkorn Academic
Advancement into Its 2nd Century Project, NSRF via the Program Management Unit
for Human Resources & Institutional Development, Research and Innovation
[grant numbers B05F650021, B37G660013] (Thailand) for providing computing
infrastructure that has contributed to the research results reported within
this paper. URL:www.e-science.in.th.

## References

  * Regge (1959) T. Regge, Nuovo Cim. 14, 951 (1959).
  * Collins (2009) P. D. B. Collins, _An Introduction to Regge Theory and High-Energy Physics_ , Cambridge Monographs on Mathematical Physics (Cambridge Univ. Press, Cambridge, UK, 2009), ISBN 978-0-521-11035-8.
  * Collins and Martin (1984) P. D. B. Collins and A. D. Martin, _HADRON INTERACTIONS_ (1984), ISBN 978-0-85274-768-1.
  * Gribov (2008) V. Gribov, _Strong Interactions of Hadrons at High Energies: Gribov Lectures on Theoretical Physics_ , Cambridge Monographs on Particle Physics, Nuclear Physics and Cosmology (Cambridge University Press, 2008).
  * Froissart (1961) M. Froissart, Phys. Rev. 123, 1053 (1961).
  * Lukaszuk and Martin (1967) L. Lukaszuk and A. Martin, Nuovo Cim. A 52, 122 (1967).
  * Chew and Frautschi (1962) G. F. Chew and S. C. Frautschi, Phys. Rev. Lett. 8, 41 (1962).
  * Donnachie et al. (2002) S. Donnachie, G. Dosch, P. Landshoff, and O. Nachtmann, _Pomeron Physics and QCD_ , Cambridge Monographs on Particle Physics, Nuclear Physics and Cosmology (Cambridge University Press, 2002).
  * Barone and Predazzi (2002) V. Barone and E. Predazzi, _High-Energy Particle Diffraction_ , vol. v.565 of _Texts and Monographs in Physics_ (Springer-Verlag, Berlin Heidelberg, 2002), ISBN 978-3-540-42107-8.
  * Lukaszuk and Nicolescu (1973) L. Lukaszuk and B. Nicolescu, Lett. Nuovo Cim. 8, 405 (1973).
  * Bartels et al. (2000) J. Bartels, L. N. Lipatov, and G. P. Vacca, Phys. Lett. B 477, 178 (2000), eprint hep-ph/9912423.
  * Ewerz et al. (2014) C. Ewerz, M. Maniatis, and O. Nachtmann, Annals Phys. 342, 31 (2014), eprint 1309.3478.
  * Ewerz et al. (2016) C. Ewerz, P. Lebiedowicz, O. Nachtmann, and A. Szczurek, Phys. Lett. B 763, 382 (2016), eprint 1606.08067.
  * Covolan et al. (1996) R. J. M. Covolan, J. Montanha, and K. A. Goulianos, Phys. Lett. B 389, 176 (1996).
  * Block and Halzen (2012) M. M. Block and F. Halzen, Phys. Rev. D 86, 051504 (2012), eprint 1208.4086.
  * Szanyi et al. (2020) I. Szanyi, L. Jenkovszky, R. Schicker, and V. Svintozelskyi, Nucl. Phys. A 998, 121728 (2020), eprint 1910.02494.
  * Jenkovszky (2020) L. Jenkovszky, Symmetry 12, 1784 (2020), eprint 2011.06432.
  * Broniowski et al. (2018) W. Broniowski, L. Jenkovszky, E. Ruiz Arriola, and I. Szanyi, Phys. Rev. D 98, 074012 (2018), eprint 1806.04756.
  * Csorgo and Szanyi (2021) T. Csorgo and I. Szanyi, Eur. Phys. J. C 81, 611 (2021), eprint 2005.14319.
  * Csörgő et al. (2021) T. Csörgő, T. Novak, R. Pasechnik, A. Ster, and I. Szanyi, Eur. Phys. J. C 81, 180 (2021), eprint 1912.11968.
  * Csörgő et al. (2019) T. Csörgő, R. Pasechnik, and A. Ster, Eur. Phys. J. C 79, 62 (2019), eprint 1807.02897.
  * Ster et al. (2015) A. Ster, L. Jenkovszky, and T. Csorgo, Phys. Rev. D 91, 074018 (2015), eprint 1501.03860.
  * Block and Cahn (1985) M. M. Block and R. N. Cahn, Rev. Mod. Phys. 57, 563 (1985).
  * Khoze et al. (2018) V. A. Khoze, A. D. Martin, and M. G. Ryskin, Phys. Lett. B 780, 352 (2018), eprint 1801.07065.
  * Halzen et al. (1993) F. Halzen, G. I. Krein, and A. A. Natale, Phys. Rev. D 47, 295 (1993).
  * Donnachie and Landshoff (1983) A. Donnachie and P. V. Landshoff, Phys. Lett. B 123, 345 (1983).
  * Ma et al. (2001) W.-X. Ma, A. W. Thomas, P.-N. Shen, and L.-J. Zhou, Commun. Theor. Phys. 36, 577 (2001).
  * Hu et al. (2002) Z.-H. Hu, L.-J. Zhou, W.-X. Ma, J. Zhang, and J.-F. Liu, Commun. Theor. Phys. 38, 65 (2002).
  * He et al. (2003) X.-R. He, L.-J. Zhou, and W.-X. Ma, Commun. Theor. Phys. 39, 78 (2003).
  * Hu et al. (2008) Z.-H. Hu, L.-J. Zhou, and W.-X. Ma, Commun. Theor. Phys. 49, 729 (2008).
  * Zhou et al. (2006) L.-J. Zhou, Z.-H. Hu, and W.-X. Ma, Commun. Theor. Phys. 45, 1069 (2006).
  * Lu et al. (2020) J. Lu, L.-J. Zhou, and Z.-J. Fang, Chin. Phys. C 44, 024105 (2020).
  * Domokos et al. (2009) S. K. Domokos, J. A. Harvey, and N. Mann, Phys. Rev. D 80, 126015 (2009), eprint 0907.1084.
  * Domokos et al. (2010) S. K. Domokos, J. A. Harvey, and N. Mann, Phys. Rev. D 82, 106007 (2010), eprint 1008.2963.
  * Avsar et al. (2010) E. Avsar, Y. Hatta, and T. Matsuo, JHEP 03, 037 (2010), eprint 0912.3806.
  * Hu et al. (2018) Z. Hu, B. Maddock, and N. Mann, JHEP 08, 093 (2018), eprint 1710.02463.
  * Xie et al. (2019) W. Xie, A. Watanabe, and M. Huang, JHEP 10, 053 (2019), eprint 1901.09564.
  * Burikham and Samart (2019) P. Burikham and D. Samart, Eur. Phys. J. C 79, 452 (2019), eprint 1902.05706.
  * Liu et al. (2022a) Z. Liu, W. Xie, and A. Watanabe (2022a), eprint 2210.11246.
  * Liu et al. (2022b) Z. Liu, W. Xie, F. Sun, S. Li, and A. Watanabe, Phys. Rev. D 106, 054025 (2022b), eprint 2202.08013.
  * Ballon-Bayona et al. (2017) A. Ballon-Bayona, R. Carcassés Quevedo, and M. S. Costa, JHEP 08, 085 (2017), eprint 1704.08280.
  * Iatrakis et al. (2016) I. Iatrakis, A. Ramamurti, and E. Shuryak, Phys. Rev. D 94, 045005 (2016), eprint 1602.05014.
  * Abazov et al. (2021) V. M. Abazov et al. (TOTEM, D0), Phys. Rev. Lett. 127, 062003 (2021), eprint 2012.03981.
  * Goldberger et al. (1960) M. L. Goldberger, M. T. Grisaru, S. W. MacDowell, and D. Y. Wong, Phys. Rev. 120, 2250 (1960).
  * Buttimore et al. (1999) N. H. Buttimore, B. Z. Kopeliovich, E. Leader, J. Soffer, and T. L. Trueman, Phys. Rev. D 59, 114010 (1999), eprint hep-ph/9901339.
  * Adamczyk et al. (2013) L. Adamczyk et al. (STAR), Phys. Lett. B 719, 62 (2013), eprint 1206.1928.
  * Leader and Trueman (2000) E. Leader and T. L. Trueman, Phys. Rev. D 61, 077504 (2000), eprint hep-ph/9908221.
  * Szymanowski and Zhou (2016) L. Szymanowski and J. Zhou, Phys. Lett. B 760, 249 (2016), eprint 1604.03207.
  * Hagiwara et al. (2020) Y. Hagiwara, Y. Hatta, R. Pasechnik, and J. Zhou, Eur. Phys. J. C 80, 427 (2020), eprint 2003.03680.
  * Magallanes et al. (2022) J. B. Magallanes, P. Sawasdipol, C. Pongkitivanichkul, and D. Samart, _Odderon as reggeized oddball spin-3 in $pp$ and $p\bar{p}$ elastic scattering_ (2022), eprint 2212.11169.
  * Goldberger et al. (1957) M. Goldberger, Y. Nambu, and R. Oehme, Ann. Phys. 2, 226 (1957), ISSN 0003-4916.
  * Jacob and Wick (1959) M. Jacob and G. Wick, Ann. Phys. 7, 404 (1959), ISSN 0003-4916.
  * Mandelstam (1958) S. Mandelstam, Phys. Rev. 112, 1344 (1958).
  * Peskin and Schroeder (1995) M. E. Peskin and D. V. Schroeder, _An Introduction to quantum field theory_ (Addison-Wesley, Reading, USA, 1995).


# Density and Affinity Dependent Social Segregation and Arbitrage Equilibrium
in a Multi-class Schelling Game

Venkat Venkatasubramanian<EMAIL_ADDRESS>Complex Resilient Intelligent
Systems Laboratory, Department of Chemical Engineering, Columbia University,
New York, NY 10027, U.S.A.  Jessica Shi  Department of Industrial Engineering
and Operations Research, Columbia University, New York, NY 10027, U.S.A.  Leo
Goldman  Department of Computer Science and Engineering, Columbia University,
New York, NY 10027, U.S.A.  Arun Sankar E.M  Complex Resilient Intelligent
Systems Laboratory, Department of Chemical Engineering, Columbia University,
New York, NY 10027, U.S.A.  Abhishek Sivaram  Complex Resilient Intelligent
Systems Laboratory, Department of Chemical Engineering, Columbia University,
New York, NY 10027, U.S.A.

###### Abstract

Contrary to the widely believed hypothesis that larger, denser cities promote
socioeconomic mixing, a recent study [1] reports the opposite behavior, i.e.
more segregation. Here, we present a game-theoretic model that predicts such a
density-dependent segregation outcome in both one- and two-class systems. The
model provides key insights into the analytical conditions that lead to such
behavior. Furthermore, the arbitrage equilibrium outcome implies the equality
of effective utilities among all agents. This could be interpreted as all
agents being equally ”happy” in their respective environments in our ideal
society. We believe that our model contributes towards a deeper mathematical
understanding of social dynamics and behavior, which is important as we strive
to develop more harmonious societies.

##  I Introduction

The increasing economic inequality and the related increase in socioeconomic
segregation are of concern in many societies. In the U.S., for example, this
segregation leads to wide variations in lifestyles and life outcomes. This
lack of mixing and social interactions across the socioeconomic spectrum
contributes to social tensions and political polarization. It is generally
believed that in large, dense, urbanized environments, the diverse
happenstance of social interactions among people leads to more mixing and less
segregation. However, in a recent study, Nilforoshan et al.[1] report that
large cities increase rather than reduce socioeconomic segregation, as they
can offer a greater choice of differentiated spaces aimed at specific
socioeconomic groups. They observe that ”The consistent result that larger,
denser cities are more segregated runs counter to the hypothesis that such
cities promote socioeconomic mixing by attracting diverse individuals and
constraining space in ways that oblige them to encounter one other. Our
results support the opposite hypothesis: big cities allow their inhabitants to
seek out people who are more like themselves.”

In this work, we present a mathematical model and its agent-based simulation
that demonstrate this density-dependent segregation outcome. Our work is
related to Schelling’s seminal game-theoretic model of social segregation[2].
It is an agent-based model in which agents move according to a specified
utility function that depends on their neighbors. He showed that although
agents may only have a mild preference, they still choose to segregate into
neighborhoods of similar individuals over time. Schelling’s model was
originally proposed to study social segregation phenomena but can also be
related to phase separation in chemistry and physics [3, 4, 5, 6, 7]. Hence,
in this paper, we use the terms phase separation and social segregation
interchangeably.

It turns out that this connection between two apparently very different
domains, sociology/economics and physics/chemistry, is not just coincidental.
There exists a deep connection between game theory (the mathematical framework
for modeling phenomena in sociology/economics) and statistical mechanics (the
mathematical framework for modeling equilibrium phenomena in
physics/chemistry) that was identified by Venkatasubramanian [8]. Inspired by
this insight, he developed a novel analytical framework, called statistical
teleodynamics, which is a synthesis of the central concepts and techniques of
statistical mechanics and population game theory [9, 8, 10]. In this paper, we
use this framework to model the dynamics of socioeconomic segregation, which
reveals interesting insights. We study Schelling-like systems with one-class
and two-class models to determine the analytical conditions under which they
would undergo socioeconomic segregation. We extend the analysis first
introduced by Venkatasubramanian et al.[10] to two-class systems and a wider
range of parameters.

##  II Mathematical Formulation: Pursuit of Maximum Utility and Arbitrage
Equilibrium

Our goal is to understand the fundamental principles and mechanisms of self-
organization of goal-driven social agents. Toward that end, we develop simple
models that offer an appropriate coarse-grained description of the system.
Unlike atoms and molecules, social agents do not behave precisely and
predictably. Therefore, we have deliberately tried to keep the models as
simple as possible and not be restricted by system-specific details and
nuances without losing key insights and relevance to empirical phenomena [10].

We also wish to stress that the spirit of our modeling is similar to that of
the ideal gas or the Ising model in statistical thermodynamics. Just as real
molecules are not point-like objects or are devoid of intermolecular
interactions, as assumed in the ideal gas model in statistical mechanics, we
make similar simplifying assumptions about the social agents in our model.
These can be relaxed to make them more realistic in subsequent refinements
like, for example, van der Walls did in thermodynamics. Ideal versions serve
as useful starting points and reference states for developing more
comprehensive models of self-organization in sociological systems.

The central theme in our theory is that goal-driven social agents constantly
pursue maximum utility by jockeying for better positions via self-
organization. So, we formulate the problem by first defining the _effective
utility_ , $h_{i}$, for an agent in state $i$. The effective utility is the
net sum of the benefits minus the costs for the agents. Social agents
constantly make benefit-cost trade-offs in their self-driven dynamical
behavior to improve their socioeconomic fitness, i.e. the effective utility.
This results in a delicate dynamic balance of the benefits of aggregation
versus the costs of overcrowding the agents. In other words, the benefits of
cooperation are balanced with the costs of competition.

Furthermore, driven by natural instincts, agents also balance two competing
strategies - exploitation and exploration. Exploitation takes advantage of the
opportunities in the immediate, local neighborhood. On the other hand,
exploration examines possibilities outside.

We believe that this combination of two main strategies, namely, the benefit-
cost trade-offs of the cooperation-competition strategy with an exploitation-
exploration strategy, is a fundamental and universal evolutionary mechanism
found in most living systems.

We motivate our model by initially considering a one-class system, which is
simpler to start with. Here, all agents belong to the same socioeconomic
category. As is generally done in Schelling games, we also model the space
where the agents operate as a large lattice $L$ of local neighborhoods or
blocks, each with $M$ sites that agents can occupy. There are $Q$ such blocks,
$QM$ sites, and a total of $N$ agents, with an average agent density of
$\rho_{0}=N/(QM)$. The state of an agent is defined by specifying the block
$i$ in which it is located, and the state of the system is defined by
specifying the number of agents, $N_{i}$, in block $i$, for all blocks
($i\in\\{1,\dots,Q\\}$). The density of the agents block $i$ is given by
$\rho_{i}=N_{i}/N$. Let block $i$ also have $V_{i}$ vacant sites, so
$V_{i}=M-N_{i}$. This approach is an extension of our recent model developed
for a Schelling-like game scenario [10].

We further formulate the problem by defining the effective utility, $h_{i}$,
for the one-class agents in block $i$, which agents try to maximize by moving
to better locations (i.e., other blocks), if possible. The effective utility
is the net sum of the benefits minus the costs and has four components. The
first is that an agent prefers to have more agents in its neighborhood, as
this aggregation improves its socioeconomic quality of life. Therefore, this
affinity benefit term, representing cooperation among agents, is proportional
to the number of agents in its neighborhood. We model this as $\alpha N_{i}$,
where $\alpha>0$ is a parameter.

However, this affinity benefit comes with a cost. As more and more agents
aggregate, this overcrowding results in a congestion cost term. As
Venkatasubramanian explains [8], the resulting net benefit (= benefit - cost)
function has an inverted U-like profile (see Fig. 1). This profile is found in
many net benefit vs. resource relationships in the real world. As one consumes
a resource, it initially leads to increasing net benefit; but after a point,
the cost of the resource goes up more quickly than the benefit, thus resulting
in decreasing net benefit. The simplest model of this is a quadratic function,
$\alpha N_{i}-\beta{N_{i}}^{2}$, with the quadratic term $-\beta N_{i}^{2}$
($\beta>0$) modeling the congestion cost.

Figure 1: Net benefit of a resource for $\alpha N_{i}-\beta{N_{i}}^{2}$ (
$\alpha=6$ , $\beta=1$)

Regarding exploration, agents derive a benefit by having a large number of
vacant sites to potentially move to in the future should such a need arise.
This is the instinct to explore other opportunities, as new vacant sites are
potentially new sources of socioeconomic benefits. We call this the option
benefit term, as agents have the option to move elsewhere if needed. Again,
following Venkatasubramanian [8, 11, 10], we model this as
$\gamma\ln(M-N_{i})$, $\gamma>0$. The logarithmic function captures the
diminishing utility of this option, a commonly used feature in economics and
game theory. As before, this benefit is also associated with a cost as a
result of competition among agents for these vacant sites. We model this
disutility of competition as $-\delta\ln N_{i}$, $\delta>0$ [9, 8, 11].

Combining these four components, we have the following effective utility
function $h_{i}$ for the agents in block $i$ as,

$\displaystyle h_{i}(N_{i})=\alpha N_{i}-\beta
N_{i}^{2}+\gamma\ln(M-N_{i})-\delta\ln N_{i}$ (1)

Intuitively, the first two terms in the equation model the benefit-cost trade-
off in the exploitation behavior, while the last two model a similar trade-off
in exploration.

We can set $\delta=1$ without any loss of generality. In addition, we set
$\gamma=1$ to gain analytical simplicity, but this can be relaxed later if
necessary. So we now have

$\displaystyle h_{i}(N_{i})=\alpha N_{i}-\beta N_{i}^{2}+\ln(M-N_{i})-\ln
N_{i}$ (2)

Rewriting this in terms of the density ($\rho_{i}$) of agents in block $i$,
$\rho_{i}=N_{i}/M$, and absorbing the constant $M$ into $\alpha$ and $\beta$,
we have

$\displaystyle
h_{i}(\bm{\rho})=\alpha\rho_{i}-\beta\rho_{i}^{2}+\ln(1-\rho_{i})-\ln\rho_{i}$
(3)

For simplicity, we define $u(\rho_{i})=\alpha\rho_{i}-\beta\rho_{i}^{2}$.
Therefore, the potential $\phi(\bm{\rho})$ becomes

$\displaystyle\phi(\bm{\rho})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\int
h_{i}(\bm{x})dx_{i}=\frac{M}{N}\sum_{i=1}^{n}\int h_{i}(\bm{\rho})d\rho_{i}$
$\displaystyle=$
$\displaystyle\frac{M}{N}\sum_{i=1}^{n}\int_{0}^{\rho_{i}}\left[u(\rho)+\ln(1-\rho)-\ln\rho\right]d\rho$

One can generalize the discrete formulation to a continuous one by replacing
$\rho_{i}$ by $\rho(r)$, where the density is a continuous function of radius
$r$ of the neighborhood as demonstrated by Sivaram and Venkatasubramanian [12]
in the self-organized flocking behavior of birds.

Now, according to the theory of potential games [13], an arbitrage equilibrium
is reached when the potential is maximized. We can determine the equilibrium
utility, $h^{*}$, by maximizing the potential (see [8]), but there exists a
simpler alternative that is more convenient for our purposes here. To analyze
the equilibrium behavior, we can take the simpler agent-based perspective and
exploit the fact that at equilibrium, all agents have the same effective
utility, i.e. $h_{i}=h^{*}$, for all $i$. In other words,

$\displaystyle\alpha\rho^{*}-\beta\rho^{*2}+\ln(1-\rho^{*})-\ln\rho^{*}=h^{*}$
(5)

We explore numerically the behavior of $h^{*}$ as a function of $\rho^{*}$ in
(5), as shown in Fig. 2 ( $\beta=0$, different $\alpha$). Below a threshold
value of $\alpha$ and $\beta$, the utility function is monotonic and has a
unique density (blue curve) for a given utility value. Above the threshold,
the utility is non-monotonic (green curve) and can have multiple density
values for the same utility. The red dotted line shows this. The orange curve
is the threshold behavior.

Figure 2: Effective Utility vs Density: $h$ vs $\rho$ for different $\alpha$.
The black points are the spinodal points
($\rho_{s1}=0.146,h_{s1}=2.934;\rho_{s2}=0.854,h_{s2}=5.066$). The red points
are the binodal points
($\rho_{b1}=0.021,h_{b1}=4.00;\rho_{b2}=0.979,h_{b2}=4.00$).

Note that whether all agents remain in a single phase of uniform density
dispersed throughout the region or separate into various groups is determined
by the slope $\partial h/\partial\rho\big{|}_{\rho^{*}}$, which is the second
derivative of $\phi$, $\partial^{2}\phi/\partial^{2}\rho\big{|}_{\rho^{*}}$.
This behavior is mathematically equivalent to spinodal decomposition in
thermodynamics, widely studied, for example, in the phase separation of alloys
and polymer blends [14, 15].

Although the mathematics of spinodal decomposition is equivalent to the
segregation of social agents, there are significant differences in the
mechanism of phase separation. In spinodal decomposition in a binary alloy
formed with atoms of two metals, say $A$ and $B$, spinodal decomposition
occurs when the change in free energy during the formation of $A-A$ bonds is
less than that of the formation of $A-B$ bonds. However, our utility
formulation allows phase separation even in systems with one-class agents.
This is because our formulation introduces a disutility due to crowding and
competition. These factors encourage agents to be part of large groups. The
same behavior is also exhibited by systems with multiclass agents.

In thermodynamics, the phase between the spinodal points (discussed in
Appendix A in more detail) is unstable as it corresponds to increasing the
free energy of the system, and hence the single phase splits into two phases
of different densities to lower the free energy. For the same reason, the
phases between the spinodal and binodal points are metastable, and the phases
at the binodal points are stable. A similar behavior happens here in
statistical teleodynamics as well (see Appendix A).

##  III Social Segregation in the one-class system: Parametric Regimes

In this section, we determine the parametric regime in which social
segregation occurs. To have social segregation in one-class systems, the
utility curve should be non-monotonic in nature; more specifically, it should
have a minimum and a maximum, as shown in Figure 2. Differentiating Eq.(3) we
get

$\displaystyle\frac{dh}{d\rho}=\alpha-2\beta\rho-\frac{1}{\rho}-\frac{1}{(1-\rho)}$
(6)

We know that the condition to have minimum or maximum is $\frac{dh}{d\rho}=0$.
To have a minimum and a maximum, there should be two distinct real solutions
for $\frac{dh}{d\rho}=0$ in the domain $0\leq\rho\leq 1$.

$\displaystyle\frac{dh}{d\rho}=\alpha-2\beta\rho-\frac{1}{\rho}-\frac{1}{(1-\rho)}=0$

$\displaystyle 2\beta\rho^{3}-(\alpha+2\beta)\rho^{2}+\alpha\rho-1=0$ (7)

This cubic equation can be solved numerically or using the Cardano formula.
One can search in the $\alpha$-$\beta$ space for values that produce two real
roots for Eq.(7) in the domain $0\leq\rho\leq 1$. Fig. 3, shows the
$\alpha$-$\beta$-$\rho_{0}$ region (shaded in yellow) within which phase
separation occurs. In Fig. 4, we show the 2-D slices of the yellow region of
spontaneous phase separation. For a given value of $\alpha$, $\beta$, and
$\rho_{0}$, they show the loci of the two densities (i.e., the low and high
density groups) of the corresponding equilibrium states of agents.

Figure 3: The parametric space within the yellow region is where segregation
is guaranteed to occur. Figure 4: This figure shows the cross-section of the
solid from Figure 3. Each line represents values of $\alpha$, showing the
range of densities where segregation would occur for values of $\beta$.

##  IV Agent-based simulations: One-class system

Agent-based simulations were performed for the parameters shown in Figure 2.
From the theory we had discussed earlier, we cannot expect phase separation to
occur for $\alpha\leq 4$ and $\beta=0$. However, phase separation can occur
for $\alpha>4$ and $\beta=0$ when the initial average density is within the
spinodal densities. Our agent-based simulations agree with this theoretical
prediction.

The equilibrium configuration of the agents for nine sets of simulations is
shown in Figure 5. The configurations shown in each row of Figure 5 are for a
constant density, and those shown in each column are for the same set of
$\alpha$ and $\beta$. As expected, for $\alpha=0$ (configurations A, D, and G)
and $\alpha=4$ (configurations B, E, and H), there is no phase separation at
any densities. For $\alpha=8$, phase separation does not occur for
$\rho_{0}=0.1$ (configuration C) because this density is outside the range of
spontaneous decomposition.

However, at higher densities (configurations F and I), phase separation occurs
because the densities are now within the spinodal region (the yellow region in
Fig. 3). In addition, note that the utility of both phases is equal (see Fig.
5 captions). Similar behavior was observed in cases where we kept $\alpha$
constant and varied $\beta$.

Figure 5:  Equilibrium configuration after five million time steps on a
$300\times 300$ grid, $M=50\times 50$ and $\beta=0$, for different $\rho_{0}$
and $\alpha$. Top row(A, B, C), middle row(D, E, F), and bottom row(G, H, I)
correspond to an average density of $\rho_{0}=0.1,0.25,0.5$, respectively.
Left column (A, D, G), middle column (B ,E, H), and right column (C ,F, I)
correspond to $\alpha=0,4,8$, respectively. For single-phase systems, the
final density is seen to be the same as $\rho_{0}$. For multiphase systems the
density of the two phases are F. 0.039 and 0.988 (utility of 3.512, 3.523,
respectively) I. 0.061 and 0.991 (utility of 3.223, 3.247 respectively)

The parameter regime in which phase separation occurs was theoretically
predicted and plotted in Figure 3. Our numerical simulations agree with the
theoretical predictions. The equilibrium agent configurations for a range of
$\alpha$-$\beta$ values are shown in Figure 6. The $x$ and $y$ coordinates of
each configuration indicate the values $\alpha$ and $\beta$ used in the
simulation.

At low values of $\alpha$, there is no phase separation (for any $\beta$, see
the uniform green region in Figure 6) because there is not much incentive to
socially aggregate. However, this changes as $\alpha$ increases, because the
benefit of aggregation increases (see Eq. 3). So, social segregation occurs
spontaneously, and the space divides itself into low-density (white squares)
and high-density (green squares) blocks.

This behavior can also be explained by the monotonicity of the $h-\rho$ curve.
At low values of $\alpha$ and high values of $\beta$, the $h-\rho$ curve is
monotonically decreasing. For any fixed $\beta$, larger values of $\alpha$
result in nonmonotonic behavior for the $h-\rho$ curve. Similarly, for any
fixed $\alpha$ ($\alpha>4$), small values of $\beta$ lead to the non-monotonic
behavior of the $h-\rho$ curve. This behavior can be observed in Figure 6. All
the simulations shown in Figure 6 are performed for N=22,500.

Figure 6: Metaplot showing agent configurations for a wide range of $\alpha$
and $\beta$ ($N=22500$ for all simulations)

##  V Segregation in a Two-class system: Mathematical Analysis

The formulation of a one-class system can readily be generalized to multiclass
systems. Here, we discuss a two-class system, as an example, by modifying the
vacancy terms in Eq. (2). Consider a class of agents identified as green
agents and another class of agents identified as red agents. Their utilities
can be defined as follows.

$\displaystyle\tilde{h}_{G,i}(N_{G,i},N_{R,i})$ $\displaystyle=$
$\displaystyle\tilde{\alpha}_{G}N_{G,i}-\tilde{\beta}_{G}N_{G,i}^{2}-\ln
N_{G,i}$ (8) $\displaystyle+\ln(M-N_{G,i}-N_{R,i})$
$\displaystyle\tilde{h}_{R,i}(N_{G,i},N_{R,i})$ $\displaystyle=$
$\displaystyle\tilde{\alpha}_{R}N_{R,i}-\tilde{\beta}_{R}N_{R,i}^{2}-\ln
N_{R,i}$ (9) $\displaystyle+\ln(M-N_{G,i}-N_{R,i})$

The utilities can be rewritten in terms of densities of the two classes,
$\rho_{G}$ and $\rho_{R}$, as

$\displaystyle{h}_{G,i}(\rho_{G,i},\rho_{R,i})$ $\displaystyle=$
$\displaystyle\alpha_{G}\rho_{G,i}-\beta_{G}\rho_{G,i}^{2}-\ln\rho_{G,i}$ (10)
$\displaystyle+\ln(1-\rho_{G,i}-\rho_{R,i})$
$\displaystyle{h}_{R,i}(\rho_{G,i},\rho_{R,i})$ $\displaystyle=$
$\displaystyle\alpha_{R}\rho_{R,i}-\beta_{R}\rho_{R,i}^{2}-\ln\rho_{R,i}$ (11)
$\displaystyle+\ln(1-\rho_{G,i}-\rho_{R,i})$

where $\rho_{G,i}=\frac{N_{G,i}}{M}$ , $\rho_{R,i}=\frac{N_{R,i}}{M}$,
$\alpha_{G}=\tilde{\alpha}_{G}M$, $\beta_{G}=\tilde{\beta}_{G}M^{2}$,
$\alpha_{R}=\tilde{\alpha}_{R}M$, and $\beta_{R}=\tilde{\beta}_{R}M^{2}$.
(Refer to the supplementary material for the derivation.)

For one-class systems, we showed that phase separation occurs when the
utility-density curve is nonmonotonic. For two-class systems, the utility of a
class depends not only on its density, but also on the density of the other
class as well. The utility of each class in a two-class system is shown in
Figure 7A. $\alpha=5$ and $\beta=0$ result in a nonconcave potential function.
The nonconcavity of the potential is the primary requirement for phase
separation. The conditions that result in non-concave potential are discussed
in Appendix B.

Figure 7: A. $h$ vs. $\rho_{G}$, $\rho_{R}$ ($\alpha=5,\beta=0$) B. Constant
utility surface, h = 2.38 intersecting $h$ vs. $\rho_{G}$, $\rho_{R}$. The
intersection points are marked in red. ($\alpha=5,\beta=0$)

In a one-class system, there were only three densities corresponding to
utility in the phase-separation region. However, in two-class systems,
infinite combinations of ($\rho_{G},\rho_{R}$) can exist for a specific
utility value. These infinite density pairs are marked in red in Figure 7B.
However, the interaction of the two classes of agents restricts the number of
possible combinations of coexisting densities. The detailed analysis is
provided in Appendix B.

##  VI Agent-based simulations: two-class system

Previously, we analyzed a one-class system and explained the qualitative
changes in equilibrium configurations with changes in $\alpha$, $\beta$, and
the initial densities. Similar analysis can be performed for the two-class
system.

Consider Figure 8 showing equilibrium configurations of agents for different
densities and $\alpha$s. For simplicity, $\alpha_{G}$, and $\alpha_{R}$ are
assumed to be equal (=$\alpha)$ and $\beta_{G}$ and $\beta_{R}$ are assumed to
be equal (=$\beta$) in all the simulations reported in Figure 8. $\beta$ is
fixed to 0 in all simulations. Each row of Figure 8 displays the equilibrium
configurations for fixed densities of the two classes;
$\rho_{G0}=\rho_{R0}=0.1$ in the upper row (A, B, C),
$\rho_{G0}=\rho_{R0}=0.25$ in the middle row (D, E, F) and
$\rho_{G0}=\rho_{R0}=0.4$ in the bottom row (G, H, I). Simulations were
performed for three values of $\alpha$ ($\alpha=0,2.023,5$). Each column in
Figure 8 shows the equilibrium configuration for a fixed $\alpha$; $\alpha=0$
in the left column (A,D,G), $\alpha=2.023$ in the middle column (B, E, H), and
$\alpha=5$ in the right column (C, F, I).

Like the one-class system, phase separation is not observed in the two-class
system for small values of $\alpha$ ($\alpha=0$ and $2.023$) (Figure 8A, B, D,
E, G, H) for any densities. For $\alpha=5$, phase separation does not occur at
$\rho_{G0}=\rho_{R0}=0.1$.

However, separation is observed at higher densities
($\rho_{G0}=\rho_{R0}=0.25$ and $\rho_{G0}=\rho_{R0}=0.4$). The nature of the
game-theoretic potential can explain this behavior. Phase separation is
observed when the potential-density surface is non-concave. In Appendix B, we
explain the conditions for the concavity of the potential. For concave
functions, the eigenvalues of the Hessian will be non-positive. When at least
one eigenvalue is positive, the potential is non-concave. The Hessian
eigenvalues are provided in Table 1 for each set of parameters reported in
Figure 8. It can be seen from Table 1 that both Hessian eigenvalues (EV-1 and
EV-2) are non-positive when $\alpha=0$ and $2.023$ and therefore the potential
is concave for configurations A, B, D, E, G and H in Figure 8.

Theoretically, it can be verified that the potential is concave in the entire
density domain when $\alpha\leq 2.023$. The potential is globally non-concave
when $\alpha>2.023$ and $\beta=0$. For $\alpha=5$ and $\beta=0$, the potential
is locally concave at the density (0.1,0.1) because, at this point, both
eigenvalues of the Hessian are non-positive. Therefore, phase separation is
not observed at this set of densities. However, one of the eigenvalues is
positive at densities (0.25,0.25) and (0.4,0.4) (EV-1 $=1>0$ for configuration
F and EV-1$=2.5>0$ for configuration I). Therefore, the potential is non-
concave, and as a consequence, phase separation is observed in configurations
F and I.

Figure 8:  Equilibrium configuration after five million time steps on a $300\times 300$ grid, $M=50\times 50$ and $\beta=0$, for different $\rho_{0}$ and $\alpha$. Top row(A, B, C), middle row(D, E, F), and bottom row(G, H, I) correspond to an average density of $(\rho_{G,0},\rho_{R,0})=(0.1,0.1),(0.25,0.25)$, and $(0.4,0.4)$. Left column (A, D, G), middle column (B, E, H), and right column (C ,F, I) correspond to $\alpha_{G}=\alpha_{R}=0,2.023,$ and $5$, respectively. For all simulations, $\beta_{G}=\beta_{R}=\beta=0$. For single-phase systems, the final density is seen to be the same as ($\rho_{G0},\rho_{R0}$). For multiphase systems, the densities of the two classes ($\rho^{*}_{G},\rho^{*}_{R}$) in the three phases are F. (0.012,0.863),(0.863, 0.012),(0.131,0.131) (with utilities (2.406, 2.384), (2.384, 2.406), (2.385, 2.385)) I. (0.011,0.887), (0.887,0.011), (0.154,0.154) (with utilities (2.297,2.271), (2.271,2.297), (2.271, 2.271)) Table 1: Eigenvalues of the Hessian and local concavity for configurations reported in Figure 8. Configuration |  $\alpha$ |  $\beta$ |  $\rho_{G0}$ |  $\rho_{R0}$ |  EV-1 |  EV-2 |  Local
---|---|---|---|---|---|---|---
in Figure 8 |  |  |  |  |  |  |  concavity
A |  0 |  0 |  0.1 |  0.1 |  -10 |  -12.5 |  Concave
B |  2.023 |  0 |  0.1 |  0.1 |  -7.977 |  -10.477 |  Concave
C |  5 |  0 |  0.1 |  0.1 |  -5. |  -7.5 |  Concave
D |  0 |  0 |  0.25 |  0.25 |  -4 |  -8 |  Concave
E |  2.023 |  0 |  0.25 |  0.25 |  -1.977 |  -5.977 |  Concave
F |  5 |  0 |  0.25 |  0.25 |  1 |  -3 |  Non-concave
G |  0 |  0 |  0.4 |  0.4 |  -2.5 |  -12.5 |  Concave
H |  2.023 |  0 |  0.4 |  0.4 |  -0.477 |  -10.477 |  Concave
I |  5 |  0 |  0.4 |  0.4 |  2.5 |  -7.5 |  Non-concave

The parameter regime where phase separation occurs was theoretically predicted
and plotted in Figure 13. Our numerical simulations agree with the theoretical
predictions. The equilibrium agent configurations for a range of
$\alpha$-$\beta$ values are shown in Figure 9. The $x$ and $y$ coordinates of
each configuration indicate the values $\alpha$ and $\beta$ used in the
simulation. At low values of $\alpha$, there is no phase separation because
there is not much incentive to get together. As $\alpha$ increases, the agents
come together. For any $\alpha$, the disutility due to crowding increases as
$\beta$ increases. Mathematically, at low values of $\alpha$ and high values
of $\beta$, the potential $\phi$ is concave (characterised by non-positive
eigenvalues). For any fixed $\beta$, larger values of $\alpha$ results in non-
concave potential. Similarly, for any fixed $\alpha$ ($\alpha>2.023$), small
values of $\beta$ leads to the non-concave potential. Therefore, phase
separation is expected at large values of $\alpha$ for a specific $\beta$ and
at small values of $\beta$ for a specific $\alpha$. This behavior is observed
in our simulations reported in Figure 9. All the simulations shown in Figure 9
performed for $N_{G}=N_{R}=22,500$.

Figure 9:  A metaplot showing agent configurations for a wide range of
$\alpha$ and $\beta$ ($\alpha_{G}=\alpha_{R}=\alpha$ and
$\beta_{G}=\beta_{R}=\beta$). For all simulations, $N_{G}=N_{R}=22500$.

##  VII Methodology-Python Simulations

We developed Schelling-like models for one-class and two-class systems in the
Python environment. The inputs to the program are the number of agents of each
type, $\alpha$s and $\beta$s. The number of agents in each class remains fixed
through the simulation and is given by $\bar{\rho_{G0}}\times Q\times M$ for
the green agents and $\bar{\rho_{R0}}\times Q\times M$ for the red agents. At
each iteration, an agent and a vacant coordinate outside of the agent’s
neighborhood are chosen at random. The agent’s utility in the current position
is compared to the agent’s utility if the agent were to hypothetically move to
the new vacant coordinate. If the agent’s utility is increased by moving to
the new location, the agent moves and the respective density and utility
values are updated. If not, the agent stays in its original position. Five
million such offers are made to the agents. For the appropriate parameter set,
the movement of the agents results in phase separation.

##  VIII Discussion

In this paper, we have presented a mathematical analysis of the decision
dynamics of socioeconomic segregation in one- and two-class systems using a
utility-driven game-theoretic model. For one-class systems, we showed both
analytically and in simulations that the non-monotonic structure of the
utility model causes social segregation. Depending on the parameters $\alpha$
and $\beta$, this usually occurs at higher densities. One can also predict and
explain this segregation from the perspective of a system trying to maximize
its game-theoretic potential. We have connected these two ideas by showing
that the non-monotonic nature of the utility function results in non-concave
potentials. The results of our agent-based simulations support these
theoretical predictions. The parametric regime of segregation was identified
using theory and confirmed by agent-based simulations.

More importantly, these predictions align with the results reported in the
study by Nilforoshan et al.[1], where they observed social segregation in
larger, denser cities. As noted, they observe that ”The consistent result that
larger, denser cities are more segregated runs counter to the hypothesis that
such cities promote socioeconomic mixing by attracting diverse individuals and
constraining space in ways that oblige them to encounter one other.” Our
theory predicts this density-dependent segregation, as we discussed above.

They further state, ”Our results support the opposite hypothesis: big cities
allow their inhabitants to seek out people who are more like themselves.” This
is also expected in our theory, as seen in Fig. 6. As $\alpha$ increases, the
system tends to segregate. The parameter $\alpha$ measures the social utility
that an agent enjoys by being with others. In larger cities, as Nilforoshan et
al.[1] study observes, there are many social options. This implies that
$\alpha$ will be higher in larger cities than in smaller ones. Therefore, as
seen in Fig. 6, it is not surprising that the population segregates into
different groups. In fact, this is expected by our model.

For the two-class system, as in the one-class system, social segregation
occurs when the potential is non-concave. The negative-definiteness of the
Hessian of the game-theoretic potential determines the concavity. Our
simulations of two-class systems agree with the theoretical predictions. We
have also determined the parametric regime of phase separation for the two-
class systems.

Finally, we wish to direct the reader’s attention to an interesting
observation. We note that the effective utility ($h^{*}$) enjoyed by the
agents in the low- and high-density phases is the same (see Fig. 5F and I), as
expected, because this is an arbitrage equilibrium. From a socioeconomic
perspective, this is an interesting result. Interpreting the effective utility
as a measure of ”happiness” or ”satisfaction”, we see that the two different
socioeconomic groups are equally ”happy” in their respective environments.
Even though they are segregated, they are both equally ”satisfied” with their
lifestyles in our ideal society. To put this a bit more colorfully, the person
enjoying a beverage with a small group of friends in a low-density small town
is just as ”happy” as his/her larger city counterpart in a fancy and crowded
restaurant in our ideal society. This harmonious outcome, despite segregation,
is not necessarily bad. However, since ”happiness” is such an elusive concept,
we wish to emphasize and caution that this ideal harmonious outcome, as
envisioned in our model, might not occur in real-world societies. Although
segregated populations enjoy the same effective utility in our ideal society,
this might not be the case in the real world due to various social, economic,
and political policies and constraints.

##  Appendix A Stability analysis of one-class systems

In Fig 2, we observe that for $\alpha=0$ (blue curve) and $\alpha=4$ (orange
curve), $\partial h/\partial\rho\leq 0$ (i.e. negative slope; recall that
$\partial h/\partial\rho$ = $\partial^{2}\phi/\partial^{2}\rho$). In such a
parameter regime, phase separation does not occur. However, for higher values
of $\alpha$, regions with $\partial h/\partial\rho>0$ (i.e., positive slope)
phase separation develop.

We better understand this from Fig. 10. The upper part of this figure shows
the potential ($\phi$) vs the density ($\rho$) curve (in green) for
$\alpha=6$, $\beta=0$. The plotted equation is

$\phi=\alpha\frac{\rho^{2}}{2}-\beta\frac{\rho^{3}}{3}-\rho\ln{\rho}-(1-\rho)\ln{(1-\rho)}$

Spinodal points are shown as black dots, where $\partial
h/\partial\rho\big{|}_{\rho^{*}}$ =
$\partial^{2}\phi/\partial^{2}\rho\big{|}_{\rho^{*}}=0$. The corresponding
spinodal points are also shown in Fig. 2 as black dots on the green curve
($\alpha=8$, $\beta=0$). Fig. 10 also shows the binodal points (in red,
connected by the common tangent line), where $\partial
h/\partial\rho\big{|}_{\rho^{*}}$ =
$\partial^{2}\phi/\partial^{2}\rho\big{|}_{\rho^{*}}<0$. The corresponding
binodal points are seen in Fig. 2 as red points connected by the red dotted
line. As we see, the two binodal points enjoy the same utility (4.00), which
is the arbitrage equilibrium.

The lower part of Fig. 10 shows the loci of binodal points (red curve) and
spinodal points (black curve) for different values of $\alpha$ ($\beta=0$). As
$\alpha$ changes, the binodal and spinodal points change, and for $\alpha>4$
($\beta=0$) they disappear. Within the spinodal region, shown in dark gray,
known as the miscibility gap in thermodynamics, a single phase of uniform
density is unstable and would split into two phases of different densities.
The reason is that the potential $\phi$ of a large group here is less than the
sum of the two potentials of the low-density group and the high-density group
at the binodal points.

We see this geometrically from the common tangent line connecting the binodal
points, which is above the single-phase green curve between the spinodal
points. Agents in such regions will be self-driven towards the high-density
binodal point to increase their utility. Therefore, $\phi$ increases and the
system splits into two groups of different densities.

Thus, for the green curve in Fig. 2, a self-organized, utility-driven, stable
phase separation occurs spontaneously at the binodal points (red dotted line)
at the arbitrage equilibrium. Although the miscibility gap is unstable, the
region immediately outside it, between the black and red curves, is
metastable. Beyond the red curve, one has a stable single phase of uniform
density - no phase separation here.

Figure 10: Game potential ($\phi$) curve and the spinodal and binodal points.
For $\alpha=8,\beta=0$, the spinodal densities are $0.146$ and $0.854$; the
binodal densities are $0.021$ and $0.979$.

In summary, for high values of $\alpha$ (e.g., green curve in Fig. 2),
combined with average densities in the miscibility gap, we observe the
spontaneous emergence of two phases, high- and low-density groups of agents,
at arbitrage equilibrium, socially driven by the self-actuated pursuit of
maximum utility by the agents.

Intuitively, in the high-density phase, agents derive so much more benefit
from the affinity term (due to the high $\alpha$) that it more than
compensates for the disutilities due to congestion and competition, thus
yielding a high effective utility. Similarly, in the low-density phase, the
benefits of reduced congestion and lower competition combined with increased
option benefit more than compensate for the loss of utility from the affinity
term.

Thus, every agent enjoys the same effective utility $h^{*}$ in one phase or
the other at equilibrium. This causes equilibrium because, as noted, there is
no more arbitrage incentive left for agents to switch neighborhoods.

As noted above, this analysis is mathematically equivalent to spinodal
decomposition in statistical thermodynamics, with an important difference. In
statistical thermodynamics, agents try to minimize their chemical potentials
and the free energy of the system. Here, in statistical teleodynamics, agents
try to maximize their utilities ($h_{i}$) and the game-theoretic potential
($\phi$). In thermodynamics, chemical potentials are equal at phase
equilibrium. In teleodynamics, the utilities are equal at the arbitrage
equilibrium. The parallel is striking, but not surprising, because, as
Venkatasubramanian has shown [8, 10], statistical teleodynamics is the
generalization of statistical thermodynamics for goal-driven agents.

##  Appendix B Mathematical analysis of two-class systems

Consider the case, $\alpha_{G}=\alpha_{R}=\alpha=5$ and
$\beta_{G}=\beta_{R}=\beta=0$. For now, assume that the equilibrium utility is
$h^{*}=2.38$ (this is the equilibrium utility observed in the agent-based
simulations). Due to the symmetry ($\alpha_{G}=\alpha_{R}$ and
$\beta_{G}=\beta_{R}$), both classes of agents are expected to exhibit
identical behavior. This implies that if there exists a green-rich phase,
referred to as Phase-I, with densities $\rho^{*}_{G,I}$ and $\rho^{*}_{R,I}$,
then there must also exist a red-rich phase, referred to as Phase-II with
densities $\rho^{*}_{G,II}$ and $\rho^{*}_{R,II}$ such that
$\rho^{*}_{G,I}=\rho^{*}_{R,II}$ and $\rho^{*}_{R,I}=\rho^{*}_{G,II}$.

In addition to these phases, there may also exist other phases where both
classes have the same densities. At equilibrium, the utilities of the two
classes in each phase should be equal owing to the symmetry of the parameters,
i.e.,

$\displaystyle
h^{*}(\rho^{*}_{G,i},\rho^{*}_{R,i})=\alpha_{G}\rho^{*}_{G,i}-\beta_{G}{\rho^{*}}^{2}_{G,i}-\ln{\rho^{*}}_{G,i}+\ln(1-{\rho^{*}}_{G,i}-{\rho^{*}}_{R,i})=\alpha_{R}\rho^{*}_{R,i}-\beta_{R}{\rho^{*}}^{2}_{R,i}-\ln{\rho^{*}}_{R,i}+\ln(1-{\rho^{*}}_{G,i}-{\rho^{*}}_{R,i})$
(12)

The constraint on the equilibrium densities emerging from Equation (12) is

$\displaystyle\alpha\rho^{*}_{G,i}-\beta{\rho^{*}}^{2}_{G,i}-\ln{\rho^{*}}_{G,i}=\alpha\rho^{*}_{R,i}-\beta{\rho^{*}}^{2}_{R,i}-\ln{\rho^{*}}_{R,i}$
(13)

Let us determine the densities satisfying Eq.(12). Refer to the plot of
$g(\rho)=\alpha\rho-\beta{\rho}^{2}-\ln{\rho}$ shown in Figure 11 A. It can be
observed that for a range of values of $g$, there exist two densities ($\rho$)
with the same value of $g$. Two such solutions are shown with brown and orange
points. The two densities represented by brown points have the same value of
$g$. Similarly, the two densities represented by the orange points have the
same value of $g$. Among the two brown points, let $\rho_{1}$ be the density
of green agents, then $\rho_{2}$ will be the density of red agents in that
phase (red-rich phase). Following the symmetry argument, there exists another
phase (green-rich phase) with density $\rho_{1}$ for the red agents and
density $\rho_{2}$ for the green agents. If any such pair of densities results
in the equilibrium utility, $h^{*}$, then that can be considered as a
potential candidate for the densities in one of the phases.

Another possibility is the existence of other phases where both classes have
equal density. An example of such a pair of densities is represented by the
black point.

Figure 11: A. $g(\rho)=\alpha\rho-\beta{\rho}^{2}-\ln{\rho}$ vs. $\rho$ ($\alpha=5,\beta=0$). The marked are solutions that satisfy Eq.(12), i.e., all combinations ($\rho_{1},\rho_{2}$), ($\rho_{2},\rho_{1}$),($\rho_{3},\rho_{4}$), ($\rho_{4},\rho_{3}$) and ($\rho_{5},\rho_{5}$) result in the same utility value, $h^{*}=2.38$. $\rho_{1}=0.012,\rho_{2}=0.863,\rho_{3}=0.055,\rho_{4}=0.495,\rho_{5}=0.132$ B. Equilibrium configuration attained in the agent-based simulations for $\alpha_{G}=\alpha_{R}=5,\beta_{G}=\beta_{R}=0,N_{G}=N_{R}=22,500$ at the end of 5 million iterations. Table 2: Derivative of the utilities of the two classes at densities reported in Figure 11. Candidate Phase |  $\rho_{G}$ |  $\rho_{R}$ |  h |  $\frac{\partial{h_{G}}}{\partial{\rho_{G}}}$ |  $\frac{\partial{h_{G}}}{\partial{\rho_{R}}}$ |  $\frac{\partial{h_{R}}}{\partial{\rho_{G}}}$ |  $\frac{\partial{h_{R}}}{\partial{\rho_{R}}}$
---|---|---|---|---|---|---|---
I |  0.012 |  0.863 |  2.38 |  -84.409 |  -7.995 |  -7.995 |  -4.154
II |  0.863 |  0.012 |  2.38 |  -4.154 |  -7.995 |  -7.995 |  -84.409
III |  0132 |  0.132 |  2.38 |  -3.951 |  -1.358 |  -1.358 |  -3.951
IV |  0.055 |  0.495 |  2.38 |  -15.512 |  -2.223 |  -2.223 |  0.759
V |  0.495 |  0.055 |  2.38 |  0.759 |  -2.223 |  -2.223 |  -15.511

All densities marked in Figure 11A are solutions to Eq.(12), i.e., the density
combinations ($\rho_{1},\rho_{2}$),
($\rho_{2},\rho_{1}$),($\rho_{3},\rho_{4}$), ($\rho_{4},\rho_{3}$) and
($\rho_{5},\rho_{5}$) result in the same utility value, $h^{*}=2.38$.
Therefore, we have five potential phases into which the mixture may separate.
Our extensive search in the density domain confirmed that there is no other
solution to Eq. (12).

Can all these five phases coexist? Any phase with $\frac{\partial
h_{G}}{\partial\rho_{G_{i}}}>0$, $\frac{\partial
h_{G}}{\partial\rho_{R_{i}}}>0$, $\frac{\partial
h_{R}}{\partial\rho_{G_{i}}}>0$ or $\frac{\partial
h_{R}}{\partial\rho_{R_{i}}}>0$ is unstable. These derivatives evaluated in
the five density pairs are provided in Table 2. Note that for Phases IV and V,
one of the utility derivatives is positive. For Phase IV,
$\frac{\partial{h_{R}}}{\partial{\rho_{R}}}>0$, indicating that the red agents
of other phases can move to Phase IV and increase their utility. These
movements will increase the density of the red agents in Phase IV, and
therefore Phase IV is unstable. Similar arguments can be applied to the green
agents in Phase V. Therefore, for a phase to coexist with other phases, the
utility derivatives must be non-positive. In summary, at equilibrium, we
anticipate phases I, II, and III to coexist.

In addition to the above requirements, considering finite number of blocks and
the definition of density of each class, the candidate density pairs must also
satisfy the constraint that the number of blocks belonging to each phase must
be an integer. Consider the separation of a mixture into three phases. Let the
densities of the two classes in the three phases be represented as
$\rho^{*}_{G,I}$,$\rho^{*}_{R,I}$,$\rho^{*}_{G,II}$,$\rho^{*}_{R,II}$,
$\rho^{*}_{G,III}$ and $\rho^{*}_{R,III}$. The corresponding number of blocks
in each phase is represented by $n_{I}$,$n_{II}$ and $n_{III}$. Then the
conservation of the number of agents requires the following constraints to be
satisfied.

$\displaystyle
M\left(n_{I}\rho^{*}_{G,I}+n_{II}\rho^{*}_{G,II}+n_{III}\rho^{*}_{G,II}\right)$
$\displaystyle=$ $\displaystyle N_{G}$ (14) $\displaystyle
M(n_{I}\rho^{*}_{R,I}+n_{II}\rho^{*}_{R,II}+n_{III}\rho^{*}_{R,II})=N_{R}$
(15) $\displaystyle n_{I}+n_{II}+n_{III}=Q$ (16)

where $n_{I},n_{II},n_{III}\in\\{0,1,2...,Q\\}$. These are soft constraints,
meaning that minor differences in the calculated densities are expected in the
simulations to ensure $n_{I}$,$n_{II}$ and $n_{III}$ are integers.

The equilibrium configuration resulting from the simulation for $\alpha=5$ and
$\beta=0$ is provided in Figure 11B. As predicted by the theory, the
simulation results in three phases. Moreover, the densities of the two classes
in each phase observed in the simulation are exactly the same as those of the
predictions.

### Game-theoretic potential: Two-class system

The utilities of the two classes provided in Eq. (10) and (11) can be used to
compute the game-theoretic potential of the system as shown in [13]. The
potential is given by

$\displaystyle\phi$
$\displaystyle=\sum_{i=1}^{Q}\Bigg{(}\alpha_{G}\frac{\rho_{G,i}^{2}}{2}-\beta_{G}\frac{\rho_{G,i}^{3}}{3}-\rho_{G,i}\
\ln\rho_{G,i}$ (17)
$\displaystyle+\alpha_{R}\frac{\rho_{R,i}^{2}}{2}-\beta_{R}\frac{\rho_{R,i}^{3}}{3}-\rho_{R,i}\
\ln\rho_{R,i}$
$\displaystyle-(1-\rho_{G,i}-\rho_{R,i})\ln(1-\rho_{G,i}-\rho_{R,i})\Bigg{)}$

### Parametric regime: Two-class system

For the one-class system, the parametric regime for phase separation was
calculated from the derivative of the utility curve. For two-class systems,
there are two independent densities, and therefore, one needs to consider the
derivatives of utility with respect to both densities. In other words, one
needs to analyze the Hessian of the potential to predict phase separation.
Phase separation does not occur when the potential curve is purely concave.
Phase separation can be expected for parameters that make the potential non-
concave.

The Hessian of the potential is defined as:

$H=\begin{bmatrix}\frac{\partial^{2}\phi}{\partial{\rho_{G}}^{2}}&\frac{\partial^{2}\phi}{\partial{\rho_{G}}\partial{\rho_{R}}}\\\
\\\\[-10.00002pt]
\frac{\partial^{2}\phi}{\partial{\rho_{R}}\partial{\rho_{G}}}&\frac{\partial^{2}\phi}{\partial{\rho_{R}}^{2}}\end{bmatrix}$
(18)

Substituting the potential expression, we get,

$H=\begin{bmatrix}\alpha_{G}-2\beta_{G}\rho_{G}-\frac{1}{\rho_{G}}-\frac{1}{1-\rho_{G}-\rho_{R}}&\frac{-1}{1-\rho_{G}-\rho_{R}}\\\
\\\
\frac{-1}{1-\rho_{G}-\rho_{R}}&\alpha_{R}-2\beta_{R}\rho_{R}-\frac{1}{\rho_{R}}-\frac{1}{1-\rho_{G}-\rho_{R}}\end{bmatrix}$
(19)

Figure 12: Shaded region indicates the region of densities of the two phases
where at least one eigenvalue of the Hessian matrix is positive
($\alpha_{G}=\alpha_{R}=\alpha=5$ and $\beta_{G}=\beta_{R}=\beta=0$)

. Phase separation occurs when the density of the initial phase is in the
shaded region.

Figure 13: Shaded region indicates $\alpha$-$\beta$ regime where phase
separation occurs. $\alpha_{G}=\alpha_{R}=\alpha$ and
$\beta_{G}=\beta_{R}=\beta$

.

The potential function is concave when the Hessian is negative semi-definite.
The easiest way to check the negative semi-definiteness of the Hessian is to
evaluate the sign of its eigenvalues. All the eigenvalues must be non-positive
for a negative semi-definite matrix. Since we are interested in a non-concave
potential function, our focus is on determining the parameter regime that
leads to positive eigenvalues for the Hessian matrix. Now, for any set of
parameters $\alpha$, $\beta$ that result in a non-concave potential, there
exists a region of densities where phase separation occurs. One of such
density regions is shown for $\alpha_{G}=\alpha_{R}=\alpha=5$ and
$\beta_{G}=\beta_{R}=\beta=0$ in Figure 12. These densities were identified by
evaluating the eigenvalues of the Hessian matrix. The points in the shaded
region in Figure 12 represent densities of the two classes that result in at
least one positive eigenvalue for the Hessian matrix.

Additionally, we performed a grid sweep in the $\alpha-\beta$ parameter
regime, searching for a non-concave potential function. The points in the
shaded region in Figure 13 represent ($\alpha,\beta$) values for which the
potential function is non-concave in the density domain. The potential
function is non-concave if the Hessian matrix has a positive eigenvalue for
any pair of densities in the entire density domain.

## References

  * Nilforoshan  _et al._ [2023] H. Nilforoshan, W. Looi, E. Pierson, B. Villanueva, N. Fishman, Y. Chen, J. Sholar, B. Redbird, D. Grusky, and J. Leskovec, Human mobility networks reveal increased segregation in large cities, Nature 624, 586 (2023).
  * Schelling [1969] T. C. Schelling, Models of segregation, The American Economic Review 59, 488 (1969).
  * 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).
  * Gauvin  _et al._ [2009] L. Gauvin, J. Vannimenus, and J.-P. Nadal, Phase diagram of a Schelling segregation model, The European Physical Journal B - Condensed Matter and Complex Systems 70, 293 (2009).
  * 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).
  * Avetisov  _et al._ [2018] V. Avetisov, A. Gorsky, S. Maslov, S. Nechaev, and O. Valba, Phase transitions in social networks inspired by the schelling model, Physical Review E 98, 032308 (2018).
  * 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), https://www.pnas.org/content/106/49/20622.full.pdf .
  * Venkatasubramanian [2017] V. Venkatasubramanian,  _How Much Inequality Is Fair?: Mathematical Principles of a Moral, Optimal, and Stable Capitalist Society_ (Columbia University Press, 2017).
  * Venkatasubramanian  _et al._ [2015] V. Venkatasubramanian, Y. Luo, and J. Sethuraman, How much inequality in income is fair?: A microeconomic game theoretic perspective, Physica A: Statistical Mechanics and its Applications 435, 120 (2015).
  * Venkatasubramanian  _et al._ [2022] V. Venkatasubramanian, A. Sivaram, and L. Das, A unified theory of emergent equilibrium phenomena in active and passive matter, Computers & Chemical Engineering 164, 107887 (2022).
  * Kanbur and Venkatasubramanian [2020] R. Kanbur and V. Venkatasubramanian, Occupational arbitrage equilibrium as an entropy maximizing solution, The European Physical Journal Special Topics 229, 1661 (2020).
  * Sivaram and Venkatasubramanian [2022] A. Sivaram and V. Venkatasubramanian, Arbitrage equilibrium, invariance, and the emergence of spontaneous order in the dynamics of birds flocking, ArXiv:2207.13743v4 (2022).
  * Sandholm [2010] W. H. Sandholm,  _Population games and evolutionary dynamics_ (MIT press, 2010).
  * Cahn [1961] J. W. Cahn, On spinodal decomposition, Acta metallurgica 9, 795 (1961).
  * Favvas and Mitropoulos [2008] E. Favvas and A. C. Mitropoulos, What is spinodal decomposition, Journal of Engineering Science and Technology Review 1, 25 (2008).