# Polarized opinion states in static networks driven by limited information
horizons

H. Pérez-Martínez<EMAIL_ADDRESS>F. Bauzá Mingueza  D. Soriano-Paños  J.
Gómez-Gardeñes  L.M. Floría

###### Abstract

The widespread emergence of opinion polarization is often attributed to the
rise of social media and the internet. These platforms can promote selective
exposure, leading to the formation of echo chambers where individuals are only
exposed to viewpoints that reinforce their existing beliefs. However,
experimental evidence shows that exposure to opposing views through cross-
cutting ties is common in both online and offline social contexts, which
frequently involve long-standing personal relationships. To account for these
facts, we have developed an opinion model that applies to static contact
structures. In this model, an agent’s influence over their neighbors depends
on the similarity of their opinions. Our findings suggest that polarization
can indeed emerge in such static structures and, driven by an increased
narrow-mindedness, even in presence of non-negligible cross-cutting ties.
Interestingly, the polarized opinion distributions generated by our model
closely resemble those obtained in surveys about highly polarized issues. This
has allowed us to categorize various issues based on their controversial
nature, sheding light on the factors that contribute to opinion polarization.

###### keywords:

Complex systems , Complex networks , Opinion formation

††journal: Chaos, Solitons and Fractals

[gotham] organization=GOTHAM lab, Institute for Biocomputation and Physics of
Complex Systems (BIFI), University of Zaragoza, postcode=50018, city=Zaragoza,
country=Spain

[matCond] organization=Department of Condensed Matter Physics, University of
Zaragoza, postcode=50009, city=Zaragoza, country=Spain

[fisTeor] organization=Department of Theoretical Physics, University of
Zaragoza, postcode=50009, city=Zaragoza, country=Spain

[gulb] organization=Institute Gulbenkian of Science (IGC), postcode=2780-156,
city=Oeiras, country=Portugal

##  1 Introduction

Opinion polarization is on the rise on modern societies [1, 2, 3, 4], and its
presence leads to significant impacts both in human behavior and in social
relationships. A clear example of such increasing influence of polarization
can be found in the political sphere. For instance, in [5] a decrease in
cross-cutting political talk has been recorded following a particularly harsh
confrontation between the governor of Wisconsin in 2011 and public-sector
employees. Moreover, different studies suggest that severe polarization leads
to legislative gridlock, democratic backsliding, political collapse [6], and
the decrease of cross-party cooperation [7]. Polarization also affects health
issues, biasing drug prescription [8] on highly politicized issues, such as
drug abuse or abortion, or compromising vaccination campaigns due to strong
in-group identification of non-vaccinated individuals. From a social
perspective, polarization sparks discrimination [9] and influences dating
behavior [10] to a degree comparable to discrimination based on race,
hampering social gatherings. For instance, in [11], the authors estimate a
loss of 34 million hours in family gatherings during Thanksgiving due to
partisan effects derived from the US national elections in 2016.

Scholars have identified social media and the pervasive use of the internet as
some of the factors contributing to this phenomenon [12, 13]. These platforms
can promote selective exposure as they allow for the formation of echo
chambers around polarized issues like gun control, abortion or vaccination
[14, 15] which can fuel radicalization processes [16]. Some recent studies
performed over web-browsing data suggest that online media consumption is
highly segregated in terms of political orientation and engagement [17], and
that social media or web-search engines exacerbates this tendency [18].

From a complex systems perspective, research has been devoted to studying echo
chambers and opinion segregation in social media, taking into account the
dynamical aspect of online social networks.In this sense, two major approaches
have been developed. The first one introduces contact networks gradually
changing over time, mimicking the behavior of social network in platforms like
Twitter. In these networks, links are created or broken either depending on
opinion similarity [19] or as a response to exposure to cross-cutting messages
[20, 21]. The other main approach uses temporal networks in which the
probability of interacting with any other agent of the system is proportional
to the opinion similarity [22, 23]. In both approaches, the major role of echo
chambers and filter bubbles in the process of opinion formation of an agent is
stressed.

In contrast, some studies suggest that the importance of echo chambers and
selective exposure in social media on the polarization processes of modern
societies might have been overstated [24, 25, 26, 27]. On the one hand, some
surveys have shown little connection between social media consumption and the
increase of polarization [28], which becomes more relevant within those groups
least likely to use the Internet [29]. On the other hand, echo chambers could
be less common than previously thought, as individuals rarely become isolated
from opposing points of view. For example, in [30] survey-based data about the
2000 general election in the United States prove that more than half of the
respondents reported cross-cutting political conversations. Moreover, recent
studies with Reddit data [31, 32] actually show a preference of political
interaction through cross-cutting ties. Likewise, offline social contexts,
such as the workplace [33] or friendship circles [34], usually feature cross-
cutting interactions, which can be transferred to social media due to the
close relationship between offline and online social interactions [35].

Stable interaction patterns, inherent to daily-life social relationships and
often featuring cross-cutting interactions [36], thus remain a key actor in
the emergence of polarization. Multiple opinion models applicable to static
contact structures have been previously developed to study the phenomenon of
opinion formation, including the Deffuant model [37] and the Hegselmann-Krause
model [38]. Both are well-known bounded confidence models, in which agents
only interact if the distance between their opinions is smaller than a certain
threshold. Both models assume a well-mixed population in their original
formulations, but their study have been also successfully extended to static
networks [39]. In these works, the dynamics leads to either consensus,
bipolarization or fragmentation of the opinion spectrum [40, 41], in which
agents usually cluster around a finite set of opinions, their number depending
on the parameters of the system. Agents also become completely isolated from
other points of view, as they lie outside their thresholds.

These models also have some drawbacks, as they consider that individuals give
the same importance to all agents regardless of the distance between their
opinions, provided that they lie inside their confidence threshold. To
alleviate this assumption, some notable extensions have been proposed: the
first one is the Relative Agreement Model [42], in which the thresholds can
change during the opinion formation process, and the influence of each
neighbor over an agent depends on the overlapping between their thresholds.
The second one is the Smooth Bounded Confidence Model [43], in which the
influence is assumed to decay following a gaussian distribution of variable
standard deviation. Regardless, both modifications remain conceptually bounded
confidence models, and qualitatively reproduce the same final opinion
distributions as the original works.

The former models fail to replicate what is observed in real world surveys
about highly polarized issues, in which people do not follow a limited set of
positions but rather display a wide range of opinions often resulting in a
bimodal distribution [44, 45]. Moreover, they do not capture cross-cutting
interactions that remain relevant in the context of daily life interaction
[30]. These limitations warrant for an increased effort to develop models to
better reproduce polarization in static networks. To fill this gap, here we
propose a new model of opinion formation over static contact structures. Our
model allows for communication between disagreeing individuals and consider
the homophily effect by which people tend to value more the opinions of like-
minded individuals. Phenomenologically, our model overcomes the former
limitations as polarized configurations with bimodal opinion distributions
naturally arises from the dynamics.

The paper is organized as follows: in section 2 we present the assumptions and
the equations of the model. In section 3 we analytically characterize the
proposed model in a fully connected graph, and in section 4 we present the
numerical results obtained by applying the model over different network
structures. In section 5 we show how the model allows retrieving opinion
distributions from experimental data and extract information about highly
polarized issues. Finally, in section 6 we summarize the main conclusions of
our work.

##  2 The model

We consider a system of $N$ interacting agents, whose acquaintances do not
change over time and are codified in an unweighted undirected adjacency matrix
$A_{ij}$. Agents’ opinions are characterized by a continuous variable
$x_{i}(t)$; the sign $\sigma(x_{i})$ represents the qualitative position of
the agent over a certain issue (for or against), and the modulus $|x_{i}|$,
her degree of radicalization.

For the opinion formation process, we choose to follow the approach in
previous works [22, 19] which assumes that agent’s opinions are influenced by
their neighbors, and change over time following three basic principles: (i)
agents lose memory of their opinions with time, (ii) they tend to align their
opinion with their neighbors, and (iii) they give more relevance to like-
minded neighbors. However, contrary to those models, we assume that
interaction patterns remain stable with time, thus implying a certain degree
of cross-cutting interactions regardless of the agents’ opinions, and in turn
consider that individuals do not pay the same attention to everyone whom they
interact with. Taking these mechanisms into account, the time evolution of the
opinion for agent $i$ can be expressed as:

$\dot{x}_{i}(t)=-x_{i}(t)+K\sum_{j=1}^{N}A_{ij}w_{ij}(t)\tanh{x_{j}(t)}\;,$
(1)

where $K$ is the social interaction strength, which determines the strength
that neighbors’ opinions have over the agent’s, and the weights $w_{ij}$
represent how valuable the view of neighbor $j$ is to agent $i$. We choose
this function as:

$w_{ij}(t)=\frac{(|x_{i}-x_{j}|+\delta)^{-\beta}}{\sum_{l=1}^{N}A_{il}(|x_{i}-x_{l}|+\delta)^{-\beta}}\;,$
(2)

where $\beta$ is the homophily parameter (the higher, the more influential
like-minded neighbors), and $\delta$ is a parameter included as a small noise,
taken to be $\delta=0.002K$ throughout the paper. Note that weights change
over time, fulfilling $\sum_{j}w_{ij}=1$, and that neighbors’ influence might
not be symmetric, so that $w_{ij}\neq w_{ji}$ in general. Moreover,
differently to the previously mentioned bounded-confidence models, our weights
choice ensures that communication between disagreeing individuals remains
relevant in the dynamics regardless of the distance in the opinion space.

The system develops on an $N$-dimensional phase-space formed by the opinions
of all agents, and its behavior is defined by a system of $N$ differential
equations whose asymptotic solutions give the opinions’ stable equilibrium
configurations. Although this system cannot be solved analytically in general,
we can study its behavior in some particular cases. First, we consider a
situation in which all agents hold the same opinion, $x_{i}=x\;\forall i$. In
this situation, the $N$-equation system is reduced to a single equation
$\dot{x}=-x+K\tanh{x}$. The expected opinion of all agents in the equilibrium
$x^{*}$ is then given by:

$x^{*}=K\tanh{x^{*}}.$ (3)

For $K\leq 1$, all agents reach $x=x^{*}=0$ and the system arrives to a
neutral consensus state. For $K>1$, all agents reach non-zero final opinions,
ending up in a radicalized state. Therefore, there exists a second order
transition in $K=1$ characterized by the order parameter $x^{*}$. Note that
$x^{*}$ also corresponds to the maximum opinion value an agent can adopt given
a $K$ value in the general scenario. Therefore, agents’ opinions $x_{i}$ are
bounded, fulfilling $x_{i}\in[-x_{\text{max}},x_{\text{max}}]$, with
$x_{\text{max}}=x^{*}$. This result does not depend on the agent’s degree due
to the normalization of the weights.

##  3 Fully-connected graph

The homogeneity assumption is incapable of predicting the existence of
polarization, defined as the coexistence of opposing views. In order to
analytically characterize the possible emergence of polarization, we consider
a fully-connected graph, so that $k_{i}=N-1\;\forall i$, and perform a
stability analysis of a polarized configuration. Let us assume that such
polarized configuration corresponds to an equilibrium state of $n$ agents
holding opinions $x^{+}>0$ and $m$ agents holding opinions $x^{-}<0$, with
$n+m=N$, $n>2$, $m>3$. Values $x^{+}$, $x^{-}$ are given by the solution the
system of equations obtained from Eq. 1:

$x^{+}=K\frac{(n-1)\tanh{x^{+}}+m\left(\frac{\delta}{x^{+}-x^{-}+\delta}\right)^{\beta}\tanh{x^{-}}}{n-1+m\left(\frac{\delta}{x^{+}-x^{-}+\delta}\right)^{\beta}}\;,$
(4a)
$x^{-}=K\frac{n\left(\frac{\delta}{x^{+}-x^{-}+\delta}\right)^{\beta}\tanh{x^{+}}+(m-1)\tanh{x^{-}}}{n\left(\frac{\delta}{x^{+}-x^{-}+\delta}\right)^{\beta}+m-1}\;.$
(4b)

This system can be solved numerically, and its behavior for different
elections of the parameters is shown in Fig. 1(a). The equilibrium values of
$x_{+}$ and $x_{-}$ depend on the relative number of individuals with positive
and negative opinions ($n$, $m$), and the values of parameters $\beta$ and
$K$. For higher values of $\beta$, agents become increasingly isolated from
opposing points of view, increasing $|x_{\pm}|$. Higher values of $K$ also
affect the equilibrium opinions through $\delta$, resulting in increased
radicalization.

We study the stability of the presented polarized configurations by adding a
small perturbation $\epsilon$ to agent’s $i$ opinion, so that
$x_{i}=x^{-}+\epsilon$. Her opinion changes following the equation:

$\displaystyle\dot{x}_{i}$ $\displaystyle=\dot{\epsilon}=-(x^{-}+\epsilon)+$
$\displaystyle+K\left[\frac{(m-1)\tanh(x^{-})+n\tanh(x^{+})\left(\frac{\epsilon+\delta}{x^{+}-x^{-}-\epsilon+\delta}\right)^{\beta}}{(m-1)+n\left(\frac{\epsilon+\delta}{x^{+}-x^{-}-\epsilon+\delta}\right)^{\beta}}\right]\;.$
(5)

Assuming $\epsilon<<\delta$, we can approximate:

$\displaystyle\left(\frac{\delta+\epsilon}{x^{+}-x^{-}-\epsilon+\delta}\right)^{\beta}\simeq$
$\displaystyle\left(\frac{\delta}{x^{+}-x^{-}+\delta}\right)^{\beta}+\beta\frac{x^{+}-x^{-}+2\delta}{(x^{+}-x^{-}+\delta)^{2}}\left(\frac{\delta}{x^{+}-x^{-}+\delta}\right)^{\beta-1}\epsilon\;.$
(6)

Furthermore, considering Equation 4b and the approximation for small
$\epsilon$:

$\frac{a\epsilon+b}{c\epsilon+d}\simeq\frac{b}{d}+\frac{ad-
bc}{d^{2}}\epsilon\;,$ (7)

we obtain:

$\dot{\epsilon}=\epsilon\left[-1+A\left(\frac{\delta}{x^{+}-x^{-}+\delta}\right)^{\beta-1}\right]\;,$
(8)

where

$A=\frac{\beta
Kn(m-1)(\tanh{x^{+}}-\tanh{x^{-}})\frac{x^{+}-x^{-}+2\delta}{(x^{+}-x^{-}+\delta)^{2}}}{\left[m-1+n\left(\frac{\delta}{x^{+}-x^{-}+\delta}\right)^{\beta}\right]^{2}}\;.$
(9)

The stability of the polarized system thus depends on the second term on the
right hand side of Eq. 8. On the one hand, the factor multiplying $A$
decreases with $\beta$, and is equal to $1$ when $\beta=1$. On the other hand,
for large $K$ values, we can approximate $x^{+}=-x^{-}\simeq K$,
$\tanh(x^{+})\simeq 1$, and, considering an evenly polarized system, we have
$m\simeq n$. If $\delta$ is small enough, we find that $A\simeq\beta$, and
thus, $A$ grows monotonically with $\beta$. Nevertheless, it can be seen in
Fig. 1(b) that the perturbation grows if $\beta>1$ ($\dot{\epsilon}>0$), and
diminishes otherwise. Therefore, this kind of configurations becomes unstable
for $\beta>1$.

Figure 1:  (a) Values of $x_{+}$ (black), $x_{-}$ (red) obtained from Eqs. 4a
and 4b. Results are shown for $K=10$ (top, solid lines) and $K=2$ (bottom,
dashed lines). (b) Values of $\dot{\epsilon}$ as a function of $\epsilon$
given by Eq. 8 for $K=10$. The perturbation grows below the transition
($\beta<1$) destabilizing the polarized states, and diminishes above the
transition ($\beta>1$). For both panels we fix $m=51$, $n=50$, and
$\delta=0.002K$. Figure 2:  (a) Average mean opinion
$\langle|\bar{x}|\rangle_{norm}$ and (b) standard deviation
$\langle\sigma\rangle_{norm}$ obtained from 100 independent realizations for
each combination of parameters $(K,\beta)$ (color code). In each realization,
initial opinions are randomly chosen on the interval
$[-x_{\text{max}},x_{\text{max}}]$. (c) Temporal evolutions of selected agents
in the representative configurations: neutral consensus (top, $K=0.5$,
$\beta=1$), radicalization (center, $K=10$, $\beta=0.1$) and polarization
(bottom, $K=10$, $\beta=1.5$). (d) Histograms corresponding to the final
configurations of panel (c). The underlying structure of contacts throughout
the figure is an Erdös-Rényi network of $N=10^{4}$ nodes with mean degree
$\langle k\rangle=10$. For the sake of clarity, only the trajectories of $1\%$
of the agents have been represented in panel (c).

##  4 Numerical results on networked populations

In order to characterize the model’s behavior in realistic environments, in
which the agents usually do not know the opinion of each of the other
components of the system, we consider static graphs with limited information
horizons, and we solve the system of $N$ coupled equations using an explicit
fourth-order Runge-Kutta method with $dt=0.1$, starting from uniformly
distributed initial states within the range
$[-x_{\text{max}},x_{\text{max}}]$. To characterize a given equilibrium
configuration, we use the absolute value of the mean opinion of the system
$|\bar{x}|=|\sum_{i}x_{i}|/N$ and its standard deviation $\sigma$, given by
$\sigma^{2}=\sum_{i}(x_{i}-\bar{x})^{2}/N$. Both order parameters allow us to
classify configurations as neutral consensus $(|\bar{x}|=0$, $\sigma=0)$,
radicalized states $(|\bar{x}|\neq 0$, $\sigma\simeq 0)$ and polarized
configurations $(|\bar{x}|\simeq 0$, $\sigma\neq 0)$. Furthermore, to classify
intermediate situations in which $|\bar{x}|\neq 0$, $\sigma\neq 0$ we consider
that a configuration in which $|\bar{x}|<\sigma$ is polarized, and radicalized
otherwise.

Figure 3:  (a) Average mean opinion $\langle|\bar{x}|\rangle_{norm}$ and (b)
standard deviation $\langle\sigma\rangle_{norm}$ obtained from 100 independent
realizations for each combination of parameters $(K,\beta)$ (color code).
Lighter colors correspond to higher values and vice versa. In each
realization, initial opinions are randomly chosen on the interval
$[-x_{\text{max}},x_{\text{max}}]$. The underlying structure of contacts
throughout the figure is an Erdös-Rényi network of $N=9824$ nodes with mean
degree $\langle k\rangle=4$.

Figs. 2(a) and (b) show the phase diagrams obtained for an Erdös-Rényi network
of $N=10^{4}$ nodes and mean degree $\langle k\rangle=10$. The process for
obtaining these results is as follows: first, we perform 100 independent
simulations for each combination of parameters $(K,\beta)$. Then, we compute
the absolute value of the mean opinion $|\bar{x}|$ and standard deviation
$\sigma$ of each of the generated opinion distributions, and average those
results obtaining $\langle|\bar{x}|\rangle$ and $\langle\sigma\rangle$. Note
that, as the range of opinions that an agent can take depends on $K$, results
are divided by $x_{\max}$ so that they lie on the range $[0,1]$. In general,
we observe the three phases predicted: for $K<1$ both measures are 0,
resulting in neutral consensus states. For $K>1$ and $\beta>1$, the mean
opinion is 0, and the standard deviation is maximum, corresponding to
polarized states. Finally, for a certain range of parameters in the region
$K>1$ and $\beta<1$, mean opinion reaches its maximum value, and the standard
deviation remains near 0. Fig. 2(c) shows some example trajectories of the
equilibrium configurations previously mentioned, taken from the marked
combination of parameters in Fig. 2(a) and (b), and Fig. 2 shows the final
opinion distributions in each of those cases.

An additional transition appears at $\beta_{c}(K)<1$. This transition splits
the radicalized phase into two: for $\beta<\beta_{c}(K)$, configurations end
up in radicalized states, while for $\beta_{c}(K)\leq\beta<1$, some
configurations end up being radicalized, while others reach metastable
polarization states. For high enough $K$ and $\beta$, behavior becomes
independent of $K$. These results prove that polarization is indeed reachable
in networked populations, and almost always present and stable for $\beta>1$.

Fig. 2 reveals the existence of a complex interplay between the homophily
mechanism and the structure of contacts in the emergence of polarization. To
explore this, in Fig. 3 we reproduce the same results as in Fig. 2 using an
Erdös-Rényi network of $N=9824$ nodes and mean degree $\langle k\rangle=4$. It
can be seen that the transition at $\beta_{c}(K)<1$ is heavily dependent on
the network degree, as is the probability of finding polarized states. In
particular, more limited information horizons result in an increased
polarization. Moreover, it can be seen in Fig. 3 that polarization can show up
when $\beta=0$ (that is, in absence of homophily) from the network structure
alone if the necessary features are present. To further investigate this
phenomenon, we repeat our analysis by fixing $K=10$ and studying the evolution
of both order parameters as we increase the average degree of the underlying
ER network. Fig. 4(a) shows that lower degrees give rise to polarization,
which quickly decays for increasing connectivity as the information horizon
grows. Therefore, if connectivity is low enough, structural communities can
breed local radicalization, and global polarization can be formed. The former
phenomenon can also be driven by other network features determining the global
communication in the network, such as the existence of structural
correlations. To illustrate it, we consider a Watts-Strogatz graph in which we
fix the degree $k=10$ and vary the correlations with the rewiring probability
$p$, and show in Fig. 4(b) that this network can breed polarization provided a
sufficiently small $p$. Taken together, these results highlight the importance
of a limited information horizon on the rise of polarization. Note that this
phenomenon cannot be characterized in models relying on temporal networks, in
which network structure was either not present or evolved as a product of the
polarization dynamics.

Figure 4:  Average mean opinion $\langle|\bar{x}|\rangle_{norm}$ (gray) and
standard deviation $\langle\sigma\rangle_{norm}$ (red) obtained from 100
independent simulations. Results are shown with 95% confidence intervals
(shadowed regions). $K=10$ and $\beta=0$ throughout the figure. (a) Results
for an Erdös-Rényi network of $N=10^{4}$ nodes and varying mean degree
$\langle k\rangle$. (b) Results for a Watts-Strogatz network of $N=10^{4}$
nodes and $k=10$, as a function of the rewiring probability $p$. The degrees
of all nodes remain fixed in the rewiring process.

Now we analyze the polarized opinions distributions generated by our model. In
Fig. 5(a), we show how these configurations are strongly shaped by the
homophily mechanism governed by the parameter $\beta$. Specifically, when
$\beta$ is low, agents’ opinions may not reach their maximum theoretical
value, but instead remain in an intermediate standpoint. This comes in stark
contrast to what can be found in previous models over static networks [39], in
which agents tend to reach a finite limited set of final opinions. The
position of the distributions’ peaks depends on $\beta$, reaching the extreme
value $x_{\text{max}}$ for $\beta\gtrsim 0.5$. By increasing $\beta$, extreme
opinions become the most common and there is a depletion of agents holding
intermediate positions until all agents reach the maximum theoretical opinion
for $\beta\gtrsim 1$, falling in line with the stability analysis performed
before.

The metastable states occurring for $\beta_{c}(K)<\beta<1$ can reach very long
lifetimes in some cases. To show this, we represent in Fig. 5(b) the fraction
of configurations considered to be polarized at $t=100$, i.e. those for which
$|\bar{x}|<\sigma$, which remain polarized at any given time $t$. There we
check that, for these $\beta$ values, polarization quickly decays in the early
stages of the evolution, but a sizeable amount of trajectories maintain
polarized configurations for a very long time. Interestingly, this fraction
increases with the size of the system, as a higher number of individuals makes
heterogeneous local environments more likely to appear.

Figure 5:  (a) Opinion histograms of metastable polarized configurations
obtained for different values of $\beta$. (b) Fraction of polarized
configurations as a function of time for $\beta=0.5$, for an ER network of
$\langle k\rangle=10$, and $10^{4}$ nodes (black line) or $10^{3}$ nodes (gray
dashed line). To perform the study we select 100 configurations considered to
be polarized at a time $t=10^{2}$ ($10^{3}$ steps), and simulate for a total
of $10^{5}$ time units, corresponding to $10^{6}$ simulation steps. A
configuration is considered to be polarized if the absolute value of its mean
opinion $|\bar{x}|$ is bigger than the standard deviation $\sigma$. The blue
vertical line marks the time chosen to establish our criterion of polarization
in the rest of the figure, which is $10^{5}$ steps corresponding to $10^{4}$
time units. (c) Fraction of polarized configurations as a function of $\beta$
for different network sizes. Results correspond to an ER network of $\langle
k\rangle=10$, comprised of $10^{4}$ nodes (black line), and $10^{3}$ nodes
(gray dashed line). (d) Fraction of polarized configurations as a function of
$\beta$ for different mean degrees. Results correspond to an ER network of
$10^{4}$ nodes and $\langle k\rangle=6$ (light gray), $\langle k\rangle=10$
(gray), and $\langle k\rangle=14$ (black). In panels (c) and (d) 200
independent simulations are performed for each value of $\beta$. We fix $K=10$
throughout the figure.

Fixing the simulation time to $t=10^{4}$ time units, we now study how the
fraction of polarized configurations also depends on $\beta$, on the
population size $N$, and the degree distribution $P(k)$. Restricting ourselves
to $\beta<1$, we see in Fig. 5(c) that bigger networks always benefit
polarization due to the larger space available for the formation of opposite,
radicalized clusters. Furthermore, Fig. 5(d) shows that lower degrees
dramatically increase the probability of reaching polarized states by limiting
the information horizons of the agents, who become increasingly isolated.
Finally, in the former figures, we can observe that, whenever the average
connectivity $\langle k\rangle$ is large enough, polarization is the most
probable around $\beta=0.5$ regardless of the network size and degree,
corresponding to the value in which, as we stated before, the opinion
histograms’ peaks reach the maximum theoretical opinion values depleting
intermediate positions.

A plausible explanation for this is that the network structure could lose its
relevance for increasing $\beta$ due to this depletion process, as each agent
would behave as if she was embedded in a fully connected network comprised of
fully radicalized agents given by her environment. Following the argument made
for fully connected networks, small perturbations of radicalized agents in
proto-polarized configurations could provoke their de-radicalization and thus,
effectively prevent the formation of polarized states. In turn, when most
agents hold intermediate positions, the assumptions of the previous study no
longer hold, agents’ environments get richer, and the network structure could
breed polarization, i.e. generate equilibrium configurations for intermediate
values of the agent’s opinions.

###  4.1 Mixing in polarized configurations

Regardless of the network structure and features, the existence of stable
contact patterns imply that agents usually become exposed to opposing views
from disagreeing neighbors, forming cross-cutting ties that hinder the rise of
echo chambers. This kind of ties is widespread on real-world societies, in
which daily-life interactions can bring divergent views close together [46,
30, 33, 34, 35]. To measure the amount of cross-cutting interactions of a
given configuration, we define the exposure $\mathcal{E}_{i}$ of agent $i$ as
the fraction of neighbors that hold opinions of opposing sign with respect to
the agents’:

$\mathcal{E}_{i}=\frac{1}{k_{i}}\sum_{j=1}^{N}A_{ij}(1-\delta_{\sigma_{i},\sigma_{j}})\;,$
(10)

where $\sigma_{i}$ indicates the sign of agent $i$’s opinion. In the case of a
randomly distributed population evenly split between both opinions, the
expected average exposure is $\bar{\mathcal{E}}=0.5$.

These links do not necessarily imply an interaction capable of greatly
disturbing an agents’ opinion because of the importance that is given to each
neighbor, determined by the weights $w_{ij}$. To measure the influence of
cross-cutting ties over agent $i$, we define her attention $\mathcal{A}_{i}$
as the sum of her cross-cutting weights:

$\mathcal{A}_{i}=\sum_{j=1}^{N}A_{ij}w_{ij}(1-\delta_{\sigma_{i},\sigma_{j}})\;.$
(11)

As shown in Fig. 6(a), the average exposure $\bar{\mathcal{E}}$ is
significantly smaller than $0.5$ when $\beta<1$, indicating an aggregation of
like-minded individuals in homophilic relationships. These opinion segregation
phenomena can bear profound consequences for the propagation of diseases [47,
48, 49] or misinformation [50, 51]. However, for $\beta>1$,
$\bar{\mathcal{E}}$ increases monotonically towards the expected result of a
well-mixed population. On the contrary, the average attention
$\bar{\mathcal{A}}$ is a monotonically decreasing function of $\beta$,
vanishing around $\beta=1$. Nevertheless, these results show that
communication between disagreeing peers is still present even in polarized
configurations, mimicking the behavior observed in real-world situations [30].

As expected, higher homophily leads to less communication between disagreeing
agents. The dependence of $\bar{\mathcal{E}}$ with $\beta$ can be therefore
explained attending to $\bar{\mathcal{A}}$: smaller values of $\beta$ sustain
information exchange in cross-cutting links, holding the agents in
intermediate positions thanks to the influence of disagreeing neighbors. This,
in turn, allows for the formation of bigger clusters of the same opinion as
originally opposing individuals maintain some degree of information exchange,
converging eventually. Higher values of $\beta$ promptly interrupt the
communication between diverging neighbors, rendering the environment’s
influence negligible and giving rise to nearly well-mixed populations. Note
that these results do not depend on the choice of $K$.

Figure 6:  (a) Average exposure (orange) and attention (purple) obtained from
100 independent simulations that reach polarization for each value of $\beta$.
(b,c,d) Exposures and attentions of the polarized configuration shown Fig.
5(a) with $\beta=0.5$. (b) Frequency of exposures and attentions. (c) Two-
dimensional histogram of individual exposure as a function of the absolute
value of the agent’s opinion. (d) Two-dimensional histogram of individual
attention as a function of the absolute value of the agent’s opinion. In all
the panels, we fix $K=10$ and consider an Erdös-Rényi network of $N=10^{4}$
nodes and $\langle k\rangle=10$.

Focusing on the local environment of the individuals, $\mathcal{E}_{i}$ span a
wide range of values regardless of the agents’ opinions, as seen in Fig. 6(c).
However, it is clear that smaller $\mathcal{E}_{i}$ are strongly linked to
intense radicalization, even though higher $\mathcal{E}_{i}$ do not guarantee
milder views. In fact, Fig. 6(d) shows that the main driver of radicalization
is the attention $\mathcal{A}_{i}$ that an agent bestows on disagreeing
neighbors: higher $\mathcal{A}_{i}$ is tantamount to weaker radicalization.
Nevertheless, even for small values of $\beta$ most agents become isolated
(albeit not completely disconnected) from opposing views, as can be seen in
Fig. 6(b). These results show that access to opposing opinions can prevent the
polarization process as long as agents remain willing to consider other points
of view [52, 16].

The linear relationship between opinion and attention observed in Fig. 6(d)
can be easily obtained from Eq. 1. Assuming $K$ high enough, the available
opinion range will be $x_{i}\in(-K,K)$, and thus, $x_{i}$ will reach higher
absolute opinion values. Then, $\tanh(x_{i})\simeq\pm 1$ for most agents, and:

$\dot{x}_{i}\simeq-
x_{i}+K\sum_{j=1}^{N}A_{ij}w_{ij}\,\text{sgn}(x_{i})=-x_{i}+K(w_{i,+}-w_{i,-})\;,$
(12)

where $w_{i,+}$ ($w_{i,-}$) represents the sum of weights with agents of
positive (negative) opinions. Recalling the normalization of weights,
$\sum_{j=1}^{N}A_{ij}w_{ij}=1$, we have that $w_{i,\pm}=1-w_{i,\mp}$, and
then, in the equilibrium:

$x_{i}=\pm K(1-2w_{i,\mp})\;,$ (13)

unveiling the linear relationship. Following this equation, it is clear that,
if $w_{i,+}>0.5$ ($w_{i,-}>0.5$), then $x_{i}>0$ ($x_{i}<0$), showing that
agents bestow a greater amount of attention to like minded neighbors. However,
this relationship only holds for $x_{j}$ high enough, so that
$\tanh(x_{j})\simeq\pm 1$. If this is not true, cross-cutting weights can
reach values greater than 0.5 (as seen in Fig. 6(d)).

##  5 Experimental Data

Finally, we apply our model to characterize the degree of polarization in the
society with respect to different topics. For this purpose, we extract the
opinions’ distributions about traditionally polarized issues from the American
National Election Study (ANES) of 2016 [45]. The ANES analyzes the voting
behavior and opinions of US society during presidential elections, by
performing nation-wide surveys about multiple topics before and after the
electoral process. We consider the topics corresponding to the following
questions:

  * 1.

2010 health care law: Do you favor, oppose, or neither favor nor oppose the
health care reform law passed in 2010? This law requires all Americans to buy
health insurance and requires health insurance companies to accept everyone.
(code: V161114x)

  * 2.

Voting as duty or choice: Different people feel differently about voting. For
some, voting is a duty - they feel they should vote in every election no
matter how they feel about the candidates and parties. For others voting is a
choice - they feel free to vote or not to vote, depending on how they feel
about the candidates and parties. For you personally, is voting mainly a duty,
mainly a choice, or neither a duty nor a choice? (code: V161151x)

  * 3.

End birthright citizenship: Some people have proposed that the U.S.
Constitution should be changed so that the children of unauthorized immigrants
do not automatically get citizenship if they are born in this country. Do you
favor, oppose, or neither favor nor oppose this proposal? (code: V161194x)

  * 4.

Wall with Mexico: Do you favor, oppose, or neither favor nor oppose building a
wall on the U.S. border with Mexico? (code: V161196x)

  * 5.

Send troops to fight ISIS: Do you favor, oppose, or neither favor nor oppose
the U.S. sending ground troops to fight Islamic militants, such as ISIS, in
Iraq and Syria? (code: V161213x)

  * 6.

Allow to refuse service to same-sex couples: Do you think business owners who
provide wedding-related services should be allowed to refuse services to same-
sex couples if same-sex marriage violates their religious beliefs, or do you
think business owners should be required to provide services regardless of a
couple’s sexual orientation? (code: V161227x)

  * 7.

Transgender bathroom use: Should transgender people – that is, people who
identify themselves as the sex or gender different from the one they were born
as – have to use the bathrooms of the gender they were born as, or should they
be allowed to use the bathrooms of their identified gender? (code: V161228x)

  * 8.

Free-trade agreements: Do you favor, oppose, or neither favor nor oppose the
U.S. making free trade agreements with other countries? (code: V162176x)

  * 9.

Government spending for healthcare: Do you favor an increase, decrease, or no
change in government spending to help people pay for health insurance when
they cannot pay for it all themselves? (code: V162193x)

  * 10.

Torture for terrorists: Do you favor, oppose, or neither favor nor oppose the
U.S. government torturing people who are suspected of being terrorists, to try
to get information? (code: V162295x)

The chosen questions correspond to different periods of surveying, before and
after the election. In particular, issues “Free-trade agreements”, “Government
spending for healthcare”, and “Torture for terrorists” correspond to the post-
election survey, and the rest, to the pre-election survey. Populations for
both surveys vary slightly, and thus, different weights must be used to
extrapolate the results to the global US population. Weights for the pre-
election survey correspond to code V160201 of the data, and for the post-
election survey, to code V160202.

Opinion histograms are normalized removing all the missing data, comprised of
answers like “don’t know”, “refused to answer”, survey errors or inexistence
of survey data. Note that “don’t know” is considered to be different from
“Neither favor nor oppose” or other similar answers equidistant from the two
extreme opinions.

###  5.1 Comparison with the model

Our aim is to find the values of $\beta$ and $K$ that generate the most
similar opinion distributions to those obtained in surveys for polarized
topics. To do that, we generate multiple polarized configurations for
different values of the parameters, and compare their opinion distributions
with the experimental ones. We use the Jensen-Shannon distance [53] in order
to measure the similarity between both distributions, and find the parameter
combination that minimize that distance.

$d_{JS}=\sqrt{\frac{D(p||m)+D(q||m)}{2}}\;,$ (14)

where $D(p||m)$ represents the Kullback-Leibler divergence:

$D(p||m)=\sum_{i}p(i)\ln{\frac{p(i)}{m(i)}}\;,$ (15)

and $m$ is the pointwise mean of $p$ and $q$. Note that $d_{JS}$ is not
defined if any value of $p$ or $q$ is zero, but no polarized distribution
generated by our model in the relevant region of the parameter space fulfills
this condition. It is also necessary that both distributions have the same
number of points or boxes to compare them (most issues are comprised of seven
data points, except for the issues “Allow to refuse service to same-sex
couples” and “Transgender bathroom use” that have six data points). Therefore,
we perform the following process over the synthetic distributions: first, for
every polarized configuration generated for a given value of $K$, we divide
the opinion value of each agent by $x_{\text{max}}(K)$, so that the final
opinions lie on the range $(-1,1)$. Then, in order to map the continuous
opinion distributions generated by the model into discrete ones, agent’s
opinions are aggregated into six or seven equally spaced boxes to match the
number of distinct choices recorded for each survey.

Figure 7:  Results for the issue “Do you think business owners who provide
wedding-related services should be allowed to refuse services to same-sex
couples if same-sex marriage violates their religious beliefs?” of the ANES
2016 (codes V161227 and V161227a, aggregated in V161227x). (a) Mean Jensen-
Shannon distance (color code) between the issue and polarized configurations
generated by our model. 100 configurations are used for each combination of
parameters $(K,\beta)$. Pairs of parameters whose probability of generating
polarized configurations is less than $5\%$ are not considered. (b) Jensen-
Shannon distances from all polarized configurations corresponding to $K=10$ to
the distribution corresponding to the aforementioned issue. A second-degree
function is fitted near the minimum to find $\beta_{\text{opt}}$ (orange
line). In purple, distances used for the fit. In gray, distances not used for
the fit. (c) Responses to the considered issue given by the ANES 2016 survey.
(d) Most similar opinion distribution obtained in our simulations. In all the
panels, we consider an Erdös-Rényi network of $N=10^{4}$ nodes and $\langle
k\rangle=10$. Figure 8:  $\beta_{\text{opt}}$ estimated for multiple polarized
issues of the ANES 2016 survey. Inset figures show some of their
distributions. For the estimation process, we fix $K=10$ and consider an
Erdös-Rényi network of $N=10^{4}$ nodes and $\langle k\rangle=10$. The errors
of the estimations are obtained applying bootstrapping techniques, but they
are too small to be visible.

The shape of our opinion distributions is not dependent on $K$, which only
determines the maximum reachable opinion. To show this, we generate 100
polarized configurations for multiple combinations of parameters $(K,\beta)$
in the range $K\in[0,20]$, $\beta\in[0,1]$ and compare them to the
experimental one corresponding to the issue “Refuse service to same sex
couples” applying the method explained above. In Fig. 7(a) we represent the
mean $d_{JS}$ for each combination $(\beta,K)$, proving that there is no
dependence with $K$ of the distances, and thus, we fix $K=10$ in the following
process to estimate $\beta_{\text{opt}}$ for all issues.

In Fig. 7(b) we show all distances between the given issue and the polarized
distributions (gray and purple points) obtained for each value of $\beta$. As
can be seen, there is a clear minimum that we can estimate by adjusting a
second-degree function (orange line) in the vicinity of the smallest distance
(purple points). The error of the adjustment is obtained applying
bootstrapping techniques, but it is too small to be visible. For the sake of
illustration, we represent in Fig. 7(c) the experimental distribution for the
selected issue, and in Fig. 7(d), the most similar synthetic distribution,
showing a great qualitative agreement among them.

Following our previous theoretical analysis of the model, $\beta_{opt}$
provides us with a way of sorting different issues by their degree of
polarization. Results for all considered issues are shown in Fig. 8. Issues
with higher $\beta_{\text{opt}}$ show more extreme polarization, as agents for
or against them tend to choose more extreme positions. Typical partisan
topics, like the construction of a wall with Mexico, the 2010 health care law
(commonly known as Obamacare), or the right of obtaining the citizenship by
birth, follow this behavior. In contrast, issues that are not usually subject
of mainstream political clash, or are perceived as extraneous by most of the
population, like free trade agreements or the use of US troops to fight ISIS
in the field, seem to arouse less exacerbating reactions. It is argued that
racial resentment and fear of losing the dominant status of white working-
class males were the main drivers of Donald Trump’s victory in the US 2016
elections rather than economic issues [54, 55, 4]. This falls in line with the
opinion about free-trade agreements and government spending on healthcare
being less polarized than other racial issues.

##  6 Conclusions

In conclusion, we have proposed a new opinion formation model applicable to
static graphs, which introduces the homophily mechanism by making use of
weights representing the attention an agent bestows on disagreeing neighbors.
We have found that polarization is indeed reachable and stable in some cases
while maintaining some amount of cross-cutting relationships, a typical
situation of daily life social interaction, and that our opinion distributions
reproduce qualitatively those of highly polarized issues observed in real-
world surveys. Furthermore, we have shown that homophily is not essential for
the rise of polarization if networks possess the appropriate features: namely
low degree or high spatial correlations. These structural features favor the
creation of limited information horizons, key in the stability of polarized
states.

The phenomenon of polarization is very complex and has been discussed from
multiple perspectives, with different definitions and measurements [56].
Following a criterion established in previous works from the complex systems
perspective [38, 22], we understand polarization as the existence of two
opposing majority opinions following a bimodal-like distribution. This
definition is also compatible with some perspectives coming from the political
science field [57, 44], which considers bimodality as a necessary condition
for polarization. The model also reproduces the phenomenon of ”group
polarization” [58] coming from the field of social psychology, as whole
populations can become increasingly radicalized towards one of the extreme
opinions given the right circumstances. In this sense, our polarized
distributions fit on the concept of ”bi-polarization” [59], as opposing groups
remain polarized in antagonistic views.

There are several limitations to our approach. For the social dynamics, our
model assumes that all individuals become engaged in the issue at hand with
similar interest and willingness to discuss, ignores any predisposition of the
agents towards a particular political orientation, and neglects other aspects
like the existence of kinship ties, reputation, or professional expertise [60,
36]. From a network perspective, for simplicity, all the results here
discussed are obtained considering synthetic models to represent the contact
structure of the population. Although we expect our results to be
qualitatively robust with respect to the network choice, extending our
analysis to real contact networks would allow disentangling the role of
different network features such as degree heterogeneities, clustering or the
existence of modular structures on the emergence and stability of polarized
states. Nevertheless, our approach constitutes a step towards the study of
opinion formation in contexts in which potentially opinion-changing
interactions are bounded by other reasons apart from simple homophily,
inherent to real-world social relationships.

## Acknowledgments

We acknowledge financial support from grant PID2020-113582GB-I00 funded by
Spanish Ministerio de Ciencia e Innovación, and from Departamento de Industria
e Innovación del Gobierno de Aragón y Fondo Social Europeo through projects
no. E30_23R (COMPHYS group) and E36_23R (FENOL group). H.P.M. and F.B.
acknowledge financial support from Gobierno de Aragón through doctoral
fellowships. F.B. acknowledges financial support from grant
PGC2018-094684-B-C22 funded by Spanish Ministerio de Ciencia e Innovación.

## References

  * [1] D. Baldassarri, P. Bearman, Dynamics of Political Polarization, American Sociological Review 72 (5) (2007) 784–811.  doi:10.1177/000312240707200507.
  * [2] D. Baldassarri, A. Gelman, Partisans without Constraint: Political Polarization and Trends in American Public Opinion, American Journal of Sociology 114 (2) (2008) 408–446.  doi:10.1086/590649.
  * [3] M. Gentzkow, Polarization in 2016, Toulouse Network for Information Technology Whitepaper (2016) 1–23.
  * [4] A. Abramowitz, J. McCoy, United States: Racial Resentment, Negative Partisanship, and Polarization in Trump’s America, The ANNALS of the American Academy of Political and Social Science 681 (1) (2019) 137–156.  doi:10.1177/0002716218811309.
  * [5] C. Wells, K. J. Cramer, M. W. Wagner, G. Alvarez, L. A. Friedland, D. V. Shah, L. Bode, S. Edgerly, I. Gabay, C. Franklin, When We Stop Talking Politics: The Maintenance and Closing of Conversation in Contentious Times, Journal of Communication 67 (1) (2017) 131–157.
  * [6] J. McCoy, T. Rahman, M. Somer, Polarization and the Global Crisis of Democracy: Common Patterns, Dynamics, and Pernicious Consequences for Democratic Polities, American Behavioral Scientist 62 (1) (2018) 16–42.  doi:10.1177/0002764218759576.
  * [7] S. Iyengar, M. Krupenkin, The Strengthening of Partisan Affect, Political Psychology 39 (S1) (2018) 201–218.  doi:10.1111/pops.12487.
  * [8] E. D. Hersh, M. N. Goldenberg, Democratic and Republican physicians provide different care on politicized health issues, Proceedings of the National Academy of Sciences 113 (42) (2016) 11811–11816, publisher: Proceedings of the National Academy of Sciences.  doi:10.1073/pnas.1606609113.
  * [9] S. Iyengar, S. J. Westwood, Fear and Loathing across Party Lines: New Evidence on Group Polarization, American Journal of Political Science 59 (3) (2015) 690–707.  doi:10.1111/ajps.12152.
  * [10] G. A. Huber, N. Malhotra, Political Homophily in Social Relationships: Evidence from Online Dating Behavior, The Journal of Politics 79 (1) (2017) 269–283, publisher: The University of Chicago Press.  doi:10.1086/687533.
  * [11] M. K. Chen, R. Rohla, The effect of partisanship and political advertising on close family ties, Science 360 (6392) (2018) 1020–1024.  doi:10.1126/science.aaq1433.
  * [12] H. Farrell, The Internet’s Consequences for Politics, Annual Review of Political Science 15 (1) (2012) 28.
  * [13] H. Allcott, L. Braghieri, S. Eichmeyer, M. Gentzkow, The Welfare Effects of Social Media, American Economic Review 110 (3) (2020) 629–676.  doi:10.1257/aer.20190658.
  * [14] A. L. Schmidt, F. Zollo, A. Scala, C. Betsch, W. Quattrociocchi, Polarization of the vaccination debate on Facebook, Vaccine 36 (25) (2018) 3606–3612.  doi:10.1016/j.vaccine.2018.05.040.
  * [15] M. Cinelli, G. D. F. Morales, A. Galeazzi, W. Quattrociocchi, M. Starnini, The echo chamber effect on social media, Proceedings of the National Academy of Sciences 118 (9) (Mar. 2021).  doi:10.1073/pnas.2023301118.
  * [16] R. Levy, Social Media, News Consumption, and Polarization: Evidence from a Field Experiment, American Economic Review 111 (3) (2021) 831–870.  doi:10.1257/aer.20191777.
  * [17] K. Garimella, T. Smith, R. Weiss, R. West, Political Polarization in Online News Consumption, Proceedings of the International AAAI Conference on Web and Social Media 15 (2021) 152–162.  doi:10.1609/icwsm.v15i1.18049.
  * [18] S. Flaxman, S. Goel, J. M. Rao, Filter Bubbles, Echo Chambers, and Online News Consumption, Public Opinion Quarterly 80 (S1) (2016) 298–320.  doi:10.1093/poq/nfw006.
  * [19] F. P. Santos, Y. Lelkes, S. A. Levin, Link recommendation algorithms and dynamics of polarization in online social networks, Proceedings of the National Academy of Sciences 118 (50) (Dec. 2021).  doi:10.1073/pnas.2102141118.
  * [20] K. Sasahara, W. Chen, H. Peng, G. L. Ciampaglia, A. Flammini, F. Menczer, Social influence and unfollowing accelerate the emergence of echo chambers, Journal of Computational Social Science 4 (1) (2021) 381–402.  doi:10.1007/s42001-020-00084-7.
  * [21] C. K. Tokita, A. M. Guess, C. E. Tarnita, Polarized information ecosystems can reorganize social networks via information cascades, Proceedings of the National Academy of Sciences 118 (50) (Dec. 2021).  doi:10.1073/pnas.2102147118.
  * [22] F. Baumann, P. Lorenz-Spreen, I. M. Sokolov, M. Starnini, Modeling Echo Chambers and Polarization Dynamics in Social Networks, Physical Review Letters 124 (4) (2020) 048301.  doi:10.1103/PhysRevLett.124.048301.
  * [23] F. Baumann, P. Lorenz-Spreen, I. M. Sokolov, M. Starnini, Emergence of Polarized Ideological Opinions in Multidimensional Topic Spaces, Physical Review X 11 (1) (2021) 011012.  doi:10.1103/PhysRevX.11.011012.
  * [24] E. Dubois, G. Blank, The echo chamber is overstated: the moderating effect of political interest and diverse media, Information, Communication & Society 21 (5) (2018) 729–745.
  * [25] A. Guess, B. Lyons, B. Nyhan, J. Reifler, Avoiding the echo chamber about echo chambers: Why selective exposure to like-minded political news is less prevalent than you think, Knight Foundation, 2018.
  * [26] Y. Benkler, R. Faris, H. Roberts, Network Propaganda: Manipulation, Disinformation, and Radicalization in American Politics, Oxford University Press, 2018.
  * [27] E. Zhuravskaya, M. Petrova, R. Enikolopov, Political Effects of the Internet and Social Media, Annual Review of Economics 12 (1) (2020) 415–438.  doi:10.1146/annurev-economics-081919-050239.
  * [28] P. Barberá, How Social Media Reduces Mass Political Polarization. Evidence from Germany, Spain, and the U.S., Job Market Paper, New York University 46 (2014) 46.
  * [29] L. Boxell, M. Gentzkow, J. M. Shapiro, Greater Internet use is not associated with faster growth in political polarization among US demographic groups, Proceedings of the National Academy of Sciences 114 (40) (2017) 10612–10617.  doi:10.1073/pnas.1706588114.
  * [30] R. Huckfeldt, J. M. Mendez, T. Osborn, Disagreement, Ambivalence, and Engagement: The Political Consequences of Heterogeneous Networks, Political Psychology 25 (1) (2004) 65–95.  doi:10.1111/j.1467-9221.2004.00357.x.
  * [31] G. De Francisci Morales, C. Monti, M. Starnini, No echo in the chambers of political interactions on Reddit, Scientific Reports 11 (1) (2021) 2818.  doi:10.1038/s41598-021-81531-x.
  * [32] C. Monti, J. D’Ignazi, M. Starnini, G. D. F. Morales, Evidence of Demographic rather than Ideological Segregation in News Discussion on Reddit (Feb. 2023).  doi:10.1145/3543507.3583468.
  * [33] D. C. Mutz, J. J. Mondak, The Workplace as a Context for Cross-Cutting Political Discourse, The Journal of Politics 68 (1) (2006) 140–155.  doi:10.1111/j.1468-2508.2006.00376.x.
  * [34] S. Goel, W. Mason, D. J. Watts, Real and perceived attitude agreement in social networks, Journal of Personality and Social Psychology 99 (2010) 611–621.  doi:10.1037/a0020697.
  * [35] E. Bakshy, S. Messing, L. A. Adamic, Exposure to ideologically diverse news and opinion on Facebook, Science 348 (6239) (2015) 1130–1132, publisher: American Association for the Advancement of Science.  doi:10.1126/science.aaa1160.
  * [36] S. Offer, C. S. Fischer, Difficult People: Who Is Perceived to Be Demanding in Personal Networks and Why Are They There?, American Sociological Review 83 (1) (2018) 111–142.  doi:10.1177/0003122417737951.
  * [37] G. Deffuant, D. Neau, F. Amblard, G. Weisbuch, Mixing beliefs among interacting agents, Advances in Complex Systems 03 (01n04) (2000) 87–98, publisher: World Scientific Publishing Co.  doi:10.1142/S0219525900000078.
  * [38] R. Hegselmann, U. Krause, Opinion Dynamics and Bounded Confidence: Models, Analysis and Simulation, Journal of Artificial Societies and Social Simulation 5 (3) (2002).
  * [39] X. F. Meng, R. A. Van Gorder, M. A. Porter, Opinion formation and distribution in a bounded-confidence model on various networks, Physical Review E 97 (2) (2018) 022312, publisher: American Physical Society.  doi:10.1103/PhysRevE.97.022312.
  * [40] J. Lorenz, Continuous opinion dynamics under bounded confidence: a survey, International Journal of Modern Physics C 18 (12) (2007) 1819–1838, publisher: World Scientific Publishing Co.  doi:10.1142/S0129183107011789.
  * [41] J. Lorenz, Fixed points in models of continuous opinion dynamics under bounded confidence, arXiv:0806.1587 [physics] (Jun. 2008).  doi:10.48550/arXiv.0806.1587.
  * [42] G. Deffuant, F. Amblard, G. Weisbuch, T. Faure, How can extremism prevail? a study based on the relative agreement interaction model, Journal of artificial societies and social simulation 5 (4) (2002).
  * [43] G. Deffuant, F. Amblard, G. Weisbuch, Modelling Group Opinion Shift to Extreme : the Smooth Bounded Confidence Model, arXiv:cond-mat/0410199 (Oct. 2004).  doi:10.48550/arXiv.cond-mat/0410199.
  * [44] M. P. Fiorina, S. J. Abrams, Political Polarization in the American Public, Annual Review of Political Science 11:563-588 (2008).
  * [45] The american national election studies (www.electionstudies.org), these materials are based on work supported by the National Science Foundation under grant numbers SES 1444721, 2014-2017, the University of Michigan, and Stanford University.
  * [46] D. C. Mutz, The Consequences of Cross-Cutting Networks for Political Participation, American Journal of Political Science 46 (4) (2002) 838–855, publisher: [Midwest Political Science Association, Wiley].  doi:10.2307/3088437.
  * [47] A. Carballosa, M. Mussa-Juane, A. P. Muñuzuri, Social Opinion Influence on Epidemic Scenarios (2021).
  * [48] G. Burgio, B. Steinegger, G. Rapisardi, A. Arenas, The impact of homophily on digital proximity tracing, Physical Review Research 3 (3) (2021) 033128.  arXiv:2103.00635, doi:10.1103/PhysRevResearch.3.033128.
  * [49] G. Burgio, B. Steinegger, A. Arenas, Homophily impacts the success of vaccine roll-outs, Communications Physics 5 (1) (2022) 70.  doi:10.1038/s42005-022-00849-8.
  * [50] D. Mocanu, L. Rossi, Q. Zhang, M. Karsai, W. Quattrociocchi, Collective attention in the age of (mis)information, Computers in Human Behavior 51 (2015) 1198–1204.  doi:10.1016/j.chb.2015.01.024.
  * [51] A. Bessi, M. Coletto, G. A. Davidescu, A. Scala, G. Caldarelli, W. Quattrociocchi, Science vs Conspiracy: Collective Narratives in the Age of Misinformation, PLOS ONE 10 (2) (2015) e0118093.  doi:10.1371/journal.pone.0118093.
  * [52] C. A. Bail, L. P. Argyle, T. W. Brown, J. P. Bumpus, H. Chen, M. B. F. Hunzaker, J. Lee, M. Mann, F. Merhout, A. Volfovsky, Exposure to opposing views on social media can increase political polarization, Proceedings of the National Academy of Sciences 115 (37) (2018) 9216–9221.  doi:10.1073/pnas.1804840115.
  * [53] J. Lin, Divergence measures based on the Shannon entropy, IEEE Transactions on Information Theory 37 (1) (1991) 145–151.  doi:10.1109/18.61115.
  * [54] S. L. Morgan, J. Lee, The White Working Class and Voter Turnout in U.S. Presidential Elections, 2004 to 2016, Sociological Science 4 (2017) 656–685.  doi:10.15195/v4.a27.
  * [55] J. Sides, M. Tesler, L. Vavreck, The 2016 U.S. Election: How Trump Lost and Won, Journal of Democracy 28 (2) (2017) 34–44.  doi:10.1353/jod.2017.0022.
  * [56] Y. Lelkes, Mass Polarization: Manifestations and Measurements, Public Opinion Quarterly 80 (S1) (2016) 392–410.  doi:10.1093/poq/nfw005.
URL https://doi.org/10.1093/poq/nfw005

  * [57] P. DiMaggio, J. Evans, B. Bryson, Have American’s Social Attitudes Become More Polarized?, American Journal of Sociology 102 (3) (1996) 690–755.  doi:10.1086/230995.
  * [58] N. E. Friedkin, Choice Shift and Group Polarization, American Sociological Review 64 (6) (1999) 856–875.  doi:10.2307/2657407.
  * [59] M. Mäs, A. Flache, Differentiation without Distancing. Explaining Bi-Polarization of Opinions without Negative Influence, PLOS ONE 8 (11) (2013) e74516.  doi:10.1371/journal.pone.0074516.
  * [60] M. L. Small, Weak ties and the core discussion network: Why people regularly discuss important matters with unimportant alters, Social Networks 35 (3) (2013) 470–483.  doi:10.1016/j.socnet.2013.05.004.


Supplement: Simulating Cosmological Evolution by Quantum Quench of an Atomic
BEC

In this supplement we present additional analytical details and simulation
results on the evolutionary dynamics of quenched Bose-Einstein condensates
(BEC). We analyze each stage sequentially in Secs. I-III, emphasizing
analogies with cosmological evolution. Also, we discuss some analytical
results deduced directly from the Gross-Pitaevskii (GP) equation, thereby
providing additional insights into our numerical simulations. In the last
section, we include details of numerical parameters used for our simulations.

##  I Dynamics in the preheating stage

To begin, we address the early dynamics in the preheating stage, focusing on
an effective theory of the upper band degrees of freedom, with the effective
Hamiltonian

$\displaystyle
H_{\text{upper}}=\sum_{\bf{k}}\epsilon_{\bf{k}}\hat{a}^{\dagger}_{\bf
k}\hat{a}_{\bf
k}+\frac{g}{2V}\sum_{\bf{k},\bf{k^{\prime}},\bf{q}}{\hat{a}^{\dagger}_{\bf{k}}{\hat{a}}^{\dagger}_{\bf{k}^{\prime}}\hat{a}_{\bf{k}^{\prime}-\bf{q}}\hat{a}_{\bf{k}+\bf{q}}}.$
(S1)

Here $\epsilon_{\bf{k}}=\hbar^{2}k_{x}^{2}/2m+|J|\cos(k_{y}d)$ is the upper
band dispersion and $\hat{a}_{\mathbf{k}},\,\hat{a}_{\mathbf{k}}^{\dagger}$
are bosonic operators in this band, $J$ is the hopping parameter associated
with the periodic potential treated here in a tight-binding approximation. To
obtain $J$, we directly diagonalize the original Hamiltonian matrix based on
the implemented lattice potential with a size cutoff where a convergence is
achieved. $J$ is exponentially dependent on the lattice depth (dominated by
$V_{1}$ here) while the band gap between first and second bands are tuned by
$V_{2}$. Because of the macroscopic occupation of a pumped BEC, we focus on
the dynamics of scattering events associated with $\mathbf{k}=\mathbf{0}$
modes. This leads to a many-body interaction term,

$\displaystyle\hat{V}_{0}\equiv\frac{g}{2V}\Big{(}{\hat{a}^{\dagger}_{\mathbf{0}}}\hat{a}^{\dagger}_{\mathbf{0}}\hat{a}_{\mathbf{0}}\hat{a}_{\mathbf{0}}+\sum_{\mathbf{k}}(\hat{a}^{\dagger}_{\mathbf{k}}\hat{a}^{\dagger}_{\mathbf{-k}}\hat{a}_{\mathbf{0}}\hat{a}_{\mathbf{0}}+h.c.)+4\sum_{\mathbf{k}}\hat{a}^{\dagger}_{\mathbf{k}}\hat{a}^{\dagger}_{\mathbf{0}}\hat{a}_{\mathbf{k}}\hat{a}_{\mathbf{0}}.\Big{)}$
(S2)

It is useful to rewrite this expression using
$\hat{a}^{\dagger}_{\mathbf{0}}\hat{a}_{\mathbf{0}}=N-\sum_{\mathbf{k}\neq\mathbf{0}}\hat{a}^{\dagger}_{\mathbf{k}}\hat{a}_{\mathbf{k}}$
where the scalar $N$ is the total number of particles in the system. Using
this replacement, one finds a quadratic effective Hamiltonian describing the
scattering events between condensed and non-condensed particles Feng  _et al._
(2018),

$\displaystyle\hat{H}^{p}_{\text{eff}}=\sum_{\mathbf{k}\neq
0}\epsilon^{\prime}_{\mathbf{k}}\hat{a}^{\dagger}_{\mathbf{k}}\hat{a}_{\mathbf{k}}+\frac{gn_{0}}{2}\sum_{\mathbf{k}\neq\mathbf{k}}(\hat{a}^{\dagger}_{\mathbf{k}}\hat{a}^{\dagger}_{-\mathbf{k}}+\hat{a}_{\mathbf{k}}\hat{a}_{-\mathbf{k}}).$
(S3)

Here $\epsilon^{\prime}_{\mathbf{k}}=\epsilon_{\mathbf{k}}-J+{gn_{0}}$ ,
$n_{0}=N/V$, and we have ignored all quantum fluctuations from non-condensed
particles. This yields Eq. 3 in the main-text. In the Heisenberg picture, we
address the time evolution of operators under this effective Hamiltonian,

$\displaystyle\hat{a}_{\mathbf{k}}(t)=e^{i\hat{H}^{p}_{\text{eff}}t{/\hbar}}\hat{a}_{\mathbf{k}}e^{-i\hat{H}^{p}_{\text{eff}}t{/\hbar}}$
(S4)

From the commutation relations it is seen that

$\displaystyle\frac{d}{dt}\begin{bmatrix}{\hat{a}_{\mathbf{k}}}(t)\\\
{\hat{a}^{\dagger}_{-\mathbf{k}}}(t)\end{bmatrix}={\frac{-i}{\hbar}}\begin{bmatrix}\epsilon^{\prime}_{\mathbf{k}}&{gn_{0}}\\\
{{-}gn_{0}}&-\epsilon^{\prime}_{\bf{k}}\end{bmatrix}\begin{bmatrix}{\hat{a}_{\mathbf{k}}}(t)\\\
{\hat{a}^{\dagger}_{-\bf{k}}}(t)\end{bmatrix}.$ (S5)

Note that the $2\times 2$ matrix here is non-hermitian so that the dynamics
are non-unitary.

From the matrix eigenvalues,

$\displaystyle\omega_{\mathbf{k}}=\sqrt{\epsilon^{\prime
2}_{\mathbf{k}}-{g^{2}n^{2}_{0}}}{/\hbar}.$ (S6)

it follows that $\omega_{\mathbf{k}}$ becomes imaginary when $\epsilon^{\prime
2}_{\mathbf{k}}-{g^{2}n^{2}_{0}}<0$ which is to be associated with exponential
growth in the number of particles. The condition
$\epsilon^{\prime}_{k}=gn_{0}$ provides the boundary separating the dynamics
associated with that of conventional Bogoliubov quasi-particles and that
associated with a parametric instability leading to the growth of a ring
condensate, as discussed in the main text.

The ring center is determined by

$\hbar^{2}k_{x}^{2}/2m+J\cos(k_{y}d)=J-gn_{0}.$ (S7)

Since $\cos(k_{y}{d})$ is periodic in momentum space, this leads to a closed
ring.

Moreover, the boundary region for exponential growth corresponds to the onset
of imaginary eigenvalues; this is associated with the condition
$|\epsilon^{\prime}_{k}|=gn_{0}$. At this boundary the modes are static,
having zero frequency.

Connection to GP simulation results— From Eq.(S7), one sees that the ring
position is directly determined by the interaction strength. This is confirmed
by numerical simulations, see Fig. S1. A more quantitative analysis is
summarized in Table S1, where a scaling of the ring position and growth
exponent is shown.

Figure S1: Ring-formation during the preheating stage, as represented by a momentum-space distribution. From left to right, the interaction energy is increased as $gn_{0}=0.0075E_{R},\,0.0151E_{R},\,0.0226E_{R}$. $n_{\bf{k}}=|\psi({\bf{k}})|^{2}$. A clear scaling of ring width with the interaction energy determined by the resonance condition is seen.  $gn_{0}/h$ (Hz) | $t\,(T)$ | $k_{y}\,(\pi/d)$ | $\Delta E\,(E_{R})$
---|---|---|---
10 | 310 | 0.28 | -0.0189
20 | 160 | 0.44 | -0.0425
30 | 110 | 0.56 | -0.0620
40 | 90 | 0.75 | -0.0892
Table S1: Summary of characteristic timescale and position for the ring
formation at different interaction strengths. $E_{R}/h=1.3$ kHz, $2\pi/T=6300$
Hz. Here $k_{y}$ is identified as the smallest $y$-momentum value in the ring,
which corresponds to $k_{x}$=0 point in the equi-energy ring in Fig. S1. The
time is chosen such that the particle density at this $k_{y}$ has grown to
around 1000 (a.u.). Here the energy difference $\Delta E=E(0,k_{y})-E_{0}$ is
measured with respect to the initial energy $E_{0}=-2.07237E_{R}$ at
${\bf{k}}=0$ in the second band. One can see a scaling of $1/t$ and $\Delta E$
as $\propto gn_{0}$. Considering the error of $\Delta t\sim 10T$, $\Delta
k_{y}\sim 0.03(\pi/d)$ in numerical data, the scaling relation is reasonably
well satisfied.

##  II Dynamics in the reheating stage

The appearance of a second condensate (through ring formation) is an important
step en route to thermalized recondensation. Crucial here is that now two
condensates are present which can lead to inter-condensate scattering events.
The ring is not a replacement of the original condensate but rather represents
a fragmentation.

###  II.1 Back-reaction: balancing between $k=0$ and ring condensate

In this subsection which is associated in the cosmological context with “back-
reaction”, we study the scattering events within ring condensates, first
establishing a mechanism to enable a co-existing BEC (corresponding to
simultaneous ring and $k=0$ condensates). Ultimately in the “re-scattering”
phase, collisions between the two lead to a de-stablization.

Once a ring BEC has formed, one has to consider scattering events involved
with ring modes. Here there are four terms to consider:

$\displaystyle\hat{V}_{r,r}$ $\displaystyle\equiv$
$\displaystyle\frac{g}{2V}\left(\sum_{\mathbf{q}_{1},\mathbf{q}_{2}\in\Lambda}\sum_{\mathbf{q}_{3},\mathbf{q}_{4}\in\Lambda}+\sum_{\mathbf{q}_{1},\mathbf{q}_{2}\in\Lambda}\sum_{\mathbf{q}_{3},\mathbf{q}_{4}\notin\Lambda}+\sum_{\mathbf{q}_{1},\mathbf{q}_{2}\notin\Lambda}\sum_{\mathbf{q}_{3},\mathbf{q}_{4}\in\Lambda}+4\sum_{\mathbf{q}_{1},\mathbf{q}_{3}\in\Lambda}\sum_{\mathbf{q}_{2},\mathbf{q}_{4}\notin\Lambda}\right)$
(S8) $\displaystyle\times$
$\displaystyle{\hat{a}_{\mathbf{q}_{4}}^{\dagger}}\hat{a}^{\dagger}_{\mathbf{q}_{3}}\hat{a}_{\mathbf{q}_{2}}\hat{a}_{\mathbf{q}_{1}}\delta(\mathbf{q}_{1}+\mathbf{q}_{2}-\mathbf{q}_{3}-\mathbf{q}_{4}).$

Since the ring is an extended object, the first summation itself contains two
types of contributions: (1) forward scattering events (such as
$\sum_{\mathbf{q}_{1},\mathbf{q}_{2}\in\Lambda}2a_{\mathbf{q}_{1}}^{\dagger}\hat{a}^{\dagger}_{\mathbf{q}_{2}}\hat{a}_{\mathbf{q}_{2}}\hat{a}_{\mathbf{q}_{1}}$);
(2) the scattering involving two incident particles with opposite momenta (
$\sum_{\mathbf{q}_{1},\mathbf{q}_{2}\in\Lambda}a_{\mathbf{q}_{1}}^{\dagger}\hat{a}^{\dagger}_{-\mathbf{q}_{1}}\hat{a}_{\mathbf{q}_{2}}\hat{a}_{-\mathbf{q}_{2}})$.
In the second/third summation, the case $\mathbf{q}_{1}=-\mathbf{q}_{2}$ is
the leading contribution while all other are sub-leading 111To observe why
$\mathbf{q}_{1}=-\mathbf{q}_{2}$ is the leading contribution, one may define
$f(Q)$ as the multiplicity for the space $q_{1}+q_{2}=Q$, i.e.,
$f(Q)=\text{Number of elements in
}\\{(q_{1},q_{2})|\,q_{1}+q_{2}=Q,\,q_{1,2}\in\Lambda\\}.$ One can check that,
for the ring-shape $\Lambda$, $f(Q)$ is on the order of $1$ if $Q\neq 0$ while
for $Q=0$, $f(Q)$ is of the order $L\gg 1$ . . In the fourth summation, the
leading contribution corresponds to $\mathbf{q}_{1}=\mathbf{q}_{3}$.

It follows that $\hat{V}_{r,r}$ can be approximated by

$\displaystyle\hat{V}_{r,r}$ $\displaystyle\simeq$
$\displaystyle\frac{g}{2V}\Bigg{(}\sum_{\mathbf{q}_{1},\mathbf{q}_{2}\in\Lambda}(2\hat{a}_{\mathbf{q}_{1}}^{\dagger}\hat{a}^{\dagger}_{\mathbf{q}_{2}}\hat{a}_{\mathbf{q}_{2}}\hat{a}_{\mathbf{q}_{1}}+\hat{a}_{\mathbf{q}_{1}}^{\dagger}\hat{a}^{\dagger}_{-\mathbf{q}_{1}}\hat{a}_{\mathbf{q}_{2}}\hat{a}_{-\mathbf{q}_{2}})$
(S9) $\displaystyle+$
$\displaystyle\sum_{\mathbf{q}\in\Lambda,\mathbf{q^{\prime}}\notin\Lambda}(\hat{a}^{\dagger}_{\mathbf{q}}\hat{a}^{\dagger}_{-\mathbf{q}}\hat{a}_{\mathbf{q^{\prime}}}\hat{a}_{-\mathbf{q^{\prime}}}+h.c.)+4\sum_{\mathbf{q}\in\Lambda,\mathbf{q^{\prime}}\notin\Lambda}\hat{a}^{\dagger}_{\mathbf{q}}\hat{a}^{\dagger}_{\mathbf{q}^{\prime}}\hat{a}_{\mathbf{q^{\prime}}}\hat{a}_{\mathbf{q}}\Bigg{)}.$

Again one can apply a mean field (M. F.) approximation and transform
${\hat{V}}_{r,r}$ to quadratic form. One may replace
$\langle{\hat{a}^{\dagger}}_{\bf{q}}\hat{a}_{\bf{q}}\rangle=N_{r}/L_{r}$ and
$\langle\hat{a}_{\bf{q}}\hat{a}_{-\bf{q}}\rangle=-iN_{r}/L_{r}$ if
${\bf{q}}\in\Lambda$. Here $L_{r}$ counts the allowed momentum modes in the
ring and $N_{r}$ is number of particles inside the ring. As discussed earlier
for the preheating stage, one thus obtains an effective Hamiltonian containing
scattering events between particles on the ring condensate, and all other
modes, among ring condensed modes and other modes,

$\displaystyle
H^{r}_{\text{eff}}=\sum_{\mathbf{k}\notin\Lambda}(\epsilon_{k}-J+gn_{0}-g{\bar{n}_{r}})\hat{a}^{\dagger}_{\mathbf{k}}\hat{a}_{\mathbf{k}}+ig{\bar{n}_{r}}(\hat{a}_{\mathbf{k}}\hat{a}_{-\mathbf{k}}-\hat{a}^{\dagger}_{\mathbf{k}}\hat{a}^{\dagger}_{-\mathbf{k}}),$
(S10)

where $\bar{n}_{r}=N_{r}/V$. The overall factor of $i$ before the quadratic
interaction term derives from the phase correlation of counter-propagating
modes in the ring. This factor may be absorbed via a redefinition of
operators, $\hat{a}_{k}\rightarrow\exp(-i\pi/4)\hat{a}_{k}$. This leads (via a
diagonalization) to a set of eigenvalues corresponding to an effective
Hamiltonian given by
$\epsilon^{r}_{k}=\sqrt{(\epsilon_{k}-J+gn_{0}-gn_{r})^{2}-g^{2}n_{r}^{2}}$.
Here we consider the limit $n_{r}=n_{0}$, where nearly all particles are
transferred to the contribution from $n_{r}$.

The key point is that the main consequence of these ring-ring scattering
events is to scatter the particles in the ring back to momenta around $k\simeq
0$. This is the counterpart of the cosmological “back-reaction”, effect which
stabilizes a co-existing BEC (corresponding to simultaneous ring and $k=0$
condensates). Ultimately in the “re-scattering” phase, collisions between the
two lead to a de-stablization. In summary, on a time scale $t\sim t_{r}$ which
corresponds to when the ring BEC has been formed, there will be a counter-flow
from the ring BEC back to modes at the band minima, until a balance of flow is
reached between two BECs.

###  II.2 Re-scattering: formation of the important cloud state

In this section, we establish a mechanism to destabilize the co-existing
BEC’s, leading to the emergence of a cloud state which represents a
proliferation of bosonic pairs.

As in the previous sections, one can write down the scattering terms between
the co-existing BECs,

$\displaystyle\hat{V}_{0,r}$ $\displaystyle\equiv$
$\displaystyle\frac{g}{V}\left(\sum_{\mathbf{k}\in\Lambda,\mathbf{q}\notin\Lambda}\hat{a}_{\mathbf{q}}^{\dagger}\hat{a}^{\dagger}_{\mathbf{k}-\mathbf{q}}\hat{a}_{\mathbf{0}}\hat{a}_{\mathbf{k}}+2\sum_{\mathbf{k}\in\Lambda}\hat{a}^{\dagger}_{\mathbf{0}}\hat{a}^{\dagger}_{\mathbf{k}}\hat{a}_{\mathbf{0}}\hat{a}_{\mathbf{k}}+2\sum_{\mathbf{k}\in\Lambda,\mathbf{q}\notin\Lambda}\hat{a}_{\mathbf{0}}^{\dagger}\hat{a}^{\dagger}_{\mathbf{q}}\hat{a}_{\mathbf{q}-\mathbf{k}}\hat{a}_{\mathbf{k}}+2\sum_{\mathbf{k}\in\Lambda,\mathbf{q}\notin\Lambda}\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}^{\dagger}_{\mathbf{q}-\mathbf{k}}\hat{a}_{\mathbf{q}}\hat{a}_{\mathbf{0}}\right)$

(a)  (b)

Figure S2:  Eigenvalue of the effective Hamiltonian in Eq. II.2. System size
is set to be 40. The figure (a) plots the imaginary part of eigenvalues. One
can observe there are a large number of eigenvalues with a small but finite
imaginary part. The figure (b) plots the real part of the eigenvalues.
Combining these two figures, one may deduce that there are a large number of
complex eigenvalues.

After a mean field approximation we arrive at an effective Hamiltonian
describing scattering events among co-existing BECs and other modes which
leads to Eq. (4) in the main text:

$\displaystyle\hat{H}^{0,r}_{\text{eff}}$ $\displaystyle\equiv$
$\displaystyle\sum_{\mathbf{k}\notin\Lambda}(\epsilon_{\mathbf{k}}-J+gn_{0}/2-g\bar{n}_{r}/2+g\bar{n}_{0}/2)\hat{a}^{\dagger}_{\mathbf{k}}\hat{a}_{\mathbf{k}}$
$\displaystyle+$ $\displaystyle
2g\sqrt{\frac{\bar{n}_{0}\bar{n}_{r}}{L_{r}}}\sum_{\mathbf{q}\in\Lambda,\mathbf{k}\notin\Lambda}(\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}^{\dagger}_{\mathbf{q}-\mathbf{k}}+h.c.)+2g\sqrt{\frac{\bar{n}_{0}\bar{n}_{r}}{L_{r}}}\sum_{\mathbf{q}\in\Lambda,\mathbf{k}\notin\Lambda}(\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}_{\mathbf{k}-\mathbf{q}}+h.c.)$

Here $\bar{n}_{0}=N_{0}/V$ is the density of $k=0$ particles. (Note that the
random phase from $\hat{a}_{\mathbf{q}}$ is unimportant as it can be absorbed
into redefined $\hat{a}_{\mathbf{k}}$ operators.) Importantly the scattering
events here are responsible for states appearing in a new momentum regime. The
Hamiltonian is intrinsically non-local in momentum space and the exact
treatment involves a diagonalization of a matrix $M$ which includes multiple
degrees of freedom. this leads to an equation of motion for $a_{k}$

$\displaystyle\frac{d}{dt}\hat{a}_{\mathbf{k}}(t)$ $\displaystyle=$
$\displaystyle{\frac{-i}{\hbar}}\Bigg{(}\epsilon^{\prime\prime}_{\mathbf{k}}\hat{a}_{\mathbf{k}}+2gn_{\text{geom}}\sum_{\mathbf{q}\in\Lambda}\frac{1}{\sqrt{L_{r}}}\left(\hat{a}_{\mathbf{q}+\mathbf{k}}(t)+\hat{a}^{\dagger}_{\mathbf{q}-\mathbf{k}}(t)\right)\Bigg{)}$
(S13)

Here
$\epsilon^{\prime\prime}_{\mathbf{k}}=\epsilon_{\mathbf{k}}-J+gn_{0}/2-g\bar{n}_{r}/2+g\bar{n}_{0}/2$
and $n_{\text{geom}}=\sqrt{\bar{n}_{0}\bar{n}_{r}}$. The closed set of
equations for all $k$ modes can be compactly written as

$\displaystyle\frac{d}{dt}\Psi(t)={\frac{-i}{\hbar}}M\Psi(t)\iff\Psi(t)=\exp\\{{-iMt/\hbar}\\}\Psi(0).$
(S14)

where the vector $\Psi$ is defined by

$\displaystyle\Psi(t)=\begin{pmatrix}\hat{a}_{\mathbf{k}_{1,1}}(t)&\hat{a}^{\dagger}_{-\mathbf{k}_{1,1}}(t)&...&...&\hat{a}_{\mathbf{k}_{L,L}}(t)&\hat{a}^{\dagger}_{-\mathbf{k}_{L,L}}(t)\end{pmatrix}^{T},\quad\mathbf{k}_{i,j}=\frac{2\pi}{L}(i,j),$
(S15)

and $L$ is the linear size of system.

Numerically diagonalizing Eq. S13, we find that the spectrum consists of a
large number of of complex eigenvalues. (See Fig S2). This can be simply
physically interpreted. In particular, since each non-local scattering term
has a small coefficient, one can deal with each pair of
$(\mathbf{q},\mathbf{k}-\mathbf{q})$ individually. One may write down the
equation of motion for a single pair of momenta (as in Eq. (5) of main text):

$\displaystyle\frac{d}{dt}\begin{bmatrix}\hat{a}_{\mathbf{k}}(t)\\\
\hat{a}^{\dagger}_{\mathbf{q}-\mathbf{k}}(t)\end{bmatrix}={\frac{-i}{\hbar}}\begin{bmatrix}\epsilon^{\prime\prime}_{\mathbf{k}}&U^{\prime}\\\
-U^{\prime}&-\epsilon^{\prime\prime}_{\mathbf{q}-\mathbf{k}}\end{bmatrix}\begin{bmatrix}\hat{a}_{\mathbf{k}}(t)\\\
\hat{a}^{\dagger}_{\mathbf{q}-\mathbf{k}}(t)\end{bmatrix}.$ (S16)

Here
$\epsilon^{\prime\prime}_{\mathbf{k}}=\epsilon_{\mathbf{k}}-J+g{n}_{0}/2-g\bar{n}_{r}/2+g\bar{n}_{0}/2$
and ${U^{\prime}}\equiv 2g\sqrt{\frac{\bar{n}_{0}\bar{n}_{r}}{L_{r}}}$. The
eigenvalue of the $2\times 2$ matrix above is given by
$\lambda=\frac{\epsilon^{\prime\prime}_{\mathbf{k}}-\epsilon^{\prime\prime}_{\mathbf{k}-\mathbf{q}}}{2}\pm\sqrt{\left(\frac{\epsilon^{\prime\prime}_{\mathbf{k}}+\epsilon^{\prime\prime}_{\mathbf{q}-\mathbf{k}}}{2}\right)^{2}-{U^{\prime}}^{2}}.$
This complex solution reflects exponential growth. Hence, the fastest growing
pair is such that
$\epsilon^{\prime\prime}_{\mathbf{k}}=-\epsilon^{\prime\prime}_{\mathbf{k}-\mathbf{q}}\iff\epsilon_{\mathbf{k}}+\epsilon_{\mathbf{k}-\mathbf{q}}=2(J-g\bar{n}_{0})$,
assuming $n_{0}=\bar{n}_{0}+\bar{n}_{r}$.

This analysis can be used to determine the shape of the cloud. Consider the
following two equations,

$\displaystyle
J\cos{(k_{y}d)}+\frac{\hbar^{2}}{2m}{k^{2}_{x}}-(J-g\bar{n}_{0})=(J-g\bar{n}_{0})-J{\cos\left[(k_{y}-q_{y})d\right]}-\frac{\hbar^{2}}{2m}{(k_{x}-q_{x})^{2}}$
$\displaystyle\text{ with }\quad
J\cos{(q_{y}d)}+\frac{\hbar^{2}}{2m}{q^{2}_{x}}=J-g{n_{0}}.$ (S17)

The first equation introduces a constraint on the kinetic energy and second
equation determines the ring shape. To further simplify these equations, one
may cancel $q_{y}$ in the first and second equations and find a single
equation below,

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
J\left(2-{\frac{\hbar^{2}q_{x}^{2}}{2mJ}-\frac{gn_{0}}{J}}\right)\cos{(k_{y}d)}-2(J-g\bar{n}_{0})$
$\displaystyle\pm$ $\displaystyle
J\sin{(k_{y}d)}\sqrt{1-{\left(1-\frac{\hbar^{2}q_{x}^{2}}{2mJ}-\frac{gn_{0}}{J}\right)^{2}}}+{\frac{\hbar^{2}}{2m}\left(2k_{x}^{2}-2k_{x}q_{x}+q_{x}^{2}\right)}\equiv
F_{{q_{x}}}(\mathbf{k}).$

It follows that the cloud shape $C$ is determined by zeros of the functions
$F_{{q_{x}}}(\mathbf{k})$, where $q_{x}$ is a free parameter. Fig. 3b in the
main text is obtained by solving for zeros of these functions numerically.

It is important to stress that after the reheating stage is completed all
remnants of ‘global’ coherence decay away at which point the cloud state is
fully formed. For the isolated quantum system we study, quantum information
cannot be lost, but it can become spread and scrambled in a more complicated
manner. Note that the cloud state is a non-thermal state. Nevertheless we
discuss next how it enables the system to thermalize.

##  III Thermalization: Boltzmann dynamics

During this last stage of equilibration, the cloud state reorganizes and the
collection of non thermal bosons contained in the cloud begin to occupy the
band minimum. This process is associated with Boltzmann dynamics, as it
involves no ‘global’ phase coherence. Since we consider weak interactions, one
may use leading order perturbative theory to characterize the processes. The
dynamical equation of each $k$-mode is given by

$\displaystyle\frac{\partial n_{\mathbf{k}}}{\partial
t}=\frac{2g^{2}}{\hbar}\int\frac{d^{2}\mathbf{k}^{\prime}d^{2}\mathbf{q}}{(2\pi)^{4}}\Gamma(\mathbf{k},\mathbf{k}^{\prime};\mathbf{k}^{\prime}+\mathbf{q},\mathbf{k}-\mathbf{q})2\pi\delta(\epsilon_{\mathbf{k}}+\epsilon_{\mathbf{k}^{\prime}}-\epsilon_{\mathbf{k}+\mathbf{q}}-\epsilon_{\mathbf{k}^{\prime}-\mathbf{q}}).$
(S19)

Here the vertex function $\Gamma$ determines how the particle number is
changed for scattering events involving the four momenta:
$\mathbf{k},\mathbf{k}^{\prime},\mathbf{k}+\mathbf{q},\mathbf{k}-\mathbf{q}^{\prime}$:

$\displaystyle\Gamma(\mathbf{k},\mathbf{k}^{\prime};\mathbf{k}+\mathbf{q},\mathbf{k}-\mathbf{q}^{\prime})=n_{\mathbf{k}-\mathbf{q}}n_{\mathbf{k}^{\prime}+\mathbf{q}}(1+n_{\mathbf{k}})(1+n_{\mathbf{k}^{\prime}})-(n_{\mathbf{k}-\mathbf{q}}+1)(n_{\mathbf{k}^{\prime}+\mathbf{q}}+1)n_{\mathbf{k}}n_{\mathbf{k}^{\prime}}.$
(S20)

The scattering event
$\mathbf{k},\mathbf{k}^{\prime}\rightarrow\mathbf{k}-\mathbf{q},\mathbf{k}^{\prime}+\mathbf{q}$
decreases $n_{\mathbf{k}}$ by unity while the reverse process similarly
increases $n_{\mathbf{k}}$. Note that
$\mathbf{k},\mathbf{k}^{\prime}\rightarrow\mathbf{k}-\mathbf{q},\mathbf{k}^{\prime}+\mathbf{q}$
and
$\mathbf{k}-\mathbf{q},\mathbf{k}^{\prime}+\mathbf{q}\rightarrow\mathbf{k},\mathbf{k}^{\prime}$
have different scattering amplitudes. Their difference leads to a cubic power
dependence on the particle number. Thus,

$\displaystyle\frac{\partial n_{\mathbf{k}}}{\partial t}=I_{\mathbf{k}}(n).$
(S21)

This equation reflects three conserved quantities: particle number, momentum,
and energy.

Additionally, this Boltzmann dynamics reflects the H-theorem which is
associated with a monotonic increase in the entropy of the Bose gas discussed
in Fig. 4 of the main text and given by

$\displaystyle S_{B}=\int
d^{2}\mathbf{k}\left[(n_{\mathbf{k}}+1)\ln(n_{\mathbf{k}}+1)-n_{\mathbf{k}}\ln
n_{\mathbf{k}}\right].$ (S22)

###  III.1 Characterization of cloud state.

To understand the equilibrating Boltzmann dynamics, it is important to provide
a description of the initial state, the cloud state $|\Psi_{\Omega}\rangle$.
This is the intermediate non-thermal state which appears after the quenching
of the condensate but before thermal equilibrium. It satisfies the following:

  * •

It is contained in a simply-connected space $\Omega$ with well defined
boundaries in momentum space. $\Omega$ is, thus, a sub-space of the two
dimensional Brillouin zone (BZ). Importantly the smaller the size of $\Omega$
the lower the kinetic energy of the system and hence the lower its effective
temperature.

  * •

For those momentum states which are present there is an effectively large
occupation number
$n_{k}(\Omega)=\langle\Psi_{\Omega}|\hat{n}_{k}|\Psi_{\Omega}\rangle$. Indeed
in our numerical simulations we find $n_{k}(\Omega)\gg 1$.

  * •

Within these occupied states there is no ‘global’ phase coherence

$\displaystyle\langle\Psi_{C}|\hat{a}_{k}\hat{a}_{-k}|\Psi_{C}\rangle\ll 1.$
(S23)

The absence of phase coherence, means that the cloud state can be treated
within a Boltzmann equation approximation.

###  III.2 The early stage of Boltzmann dynamics

The dynamics in this early stage strongly depends on these properties of the
cloud state. One may classify the 2-body collision events according to how
many particles inside and outside the cloud are involved. If there are 3 or 4
inside-the-cloud modes involved this will induce strong particle number flow,
whereas if there are only 2 such modes the flow is weaker, and with only one
or zero such modes, this contribution in the early dynamics can be ignored.

We focus, for simplicity, on a situation where the particle number
distribution in the cloud is uniform. We have tested the appropriateness of
this assumption numerically and with this assumption, one can greatly simplify
the expression for the collision integral. To this end, we define four
functions below,

$\displaystyle F_{1}(\mathbf{k})=\int_{\mathbf{k}_{2}\notin
C}\frac{d^{2}\mathbf{k}_{2}}{(2\pi)^{2}}\int_{\mathbf{k}_{3}\in
C}\frac{d^{2}\mathbf{k}_{3}}{(2\pi)^{2}}\int_{\mathbf{k}_{4}\in
C}\frac{d^{2}\mathbf{k}_{4}}{(2\pi)^{2}}2\pi\delta(\epsilon_{\mathbf{k}}+\epsilon_{\mathbf{k}_{2}}-\epsilon_{\mathbf{k}_{3}}-\epsilon_{\mathbf{k}_{4}}){\left(2\pi\right)^{2}}\delta(\mathbf{k}+\mathbf{k}_{2}-\mathbf{k}_{3}-\mathbf{k}4),$
$\displaystyle F_{2}(\mathbf{k})=\int_{\mathbf{k}_{2}\in
C}\frac{d^{2}\mathbf{k}_{2}}{(2\pi)^{2}}\int_{\mathbf{k}_{3}\notin
C}\frac{d^{2}\mathbf{k}_{3}}{(2\pi)^{2}}\int_{\mathbf{k}_{4}\in
C}\frac{d^{2}\mathbf{k}_{4}}{(2\pi)^{2}}2\pi\delta(\epsilon_{\mathbf{k}}+\epsilon_{\mathbf{k}_{2}}-\epsilon_{\mathbf{k}_{3}}-\epsilon_{\mathbf{k}_{4}}){\left(2\pi\right)^{2}}\delta(\mathbf{k}+\mathbf{k}_{2}-\mathbf{k}_{3}-\mathbf{k}4),$
$\displaystyle G(\mathbf{k})=\int_{\mathbf{k}_{2}\in
C}\frac{d^{2}\mathbf{k}_{2}}{(2\pi)^{2}}\int_{\mathbf{k}_{3}\notin
C}\frac{d^{2}\mathbf{k}_{3}}{(2\pi)^{2}}\int_{\mathbf{k}_{4}\notin
C}\frac{d^{2}\mathbf{k}_{4}}{(2\pi)^{2}}2\pi\delta(\epsilon_{\mathbf{k}}+\epsilon_{\mathbf{k}_{2}}-\epsilon_{\mathbf{k}_{3}}-\epsilon_{\mathbf{k}_{4}}){\left(2\pi\right)^{2}}\delta(\mathbf{k}+\mathbf{k}_{2}-\mathbf{k}_{3}-\mathbf{k}_{4}).$
$\displaystyle J(\mathbf{k})=\int_{\mathbf{k}_{2}\in
C}\frac{d^{2}\mathbf{k}_{2}}{(2\pi)^{2}}\int_{\mathbf{k}_{3}\in
C}\frac{d^{2}\mathbf{k}_{3}}{(2\pi)^{2}}\int_{\mathbf{k}_{4}\in
C}\frac{d^{2}\mathbf{k}_{4}}{(2\pi)^{2}}2\pi\delta(\epsilon_{\mathbf{k}}+\epsilon_{\mathbf{k}_{2}}-\epsilon_{\mathbf{k}_{3}}-\epsilon_{\mathbf{k}_{4}}){\left(2\pi\right)^{2}}\delta(\mathbf{k}+\mathbf{k}_{2}-\mathbf{k}_{3}-\mathbf{k}_{4}),$

We consider the two separate cases for $\mathbf{k}\in\Omega$ and
$\mathbf{k}\notin\Omega$. For the former the collision integral can be reduced
to be,

$\displaystyle
I_{\bf{k}}(n)\simeq\frac{2g^{2}}{\hbar}\left[F_{1}(\mathbf{k})n^{3}_{\Omega}-2F_{2}(\mathbf{k})n^{3}_{\Omega}-G(\mathbf{k})n_{\Omega}^{2}\right],$
(S24)

where only the first two terms dominate. Moreover for $\mathbf{k}\notin\Omega$
the collision integral becomes,

$I_{\mathbf{k}}(n)\simeq\frac{2g^{2}}{\hbar}\left[J(\mathbf{k})n^{3}_{\Omega}+F_{1}(\mathbf{k}){n^{2}_{\Omega}}\right].$
(S25)

Here numerically one can show that both $J(\mathbf{k})$ and
$F_{1}(\mathbf{k})$ nearly vanish except in a small region near the boundary
of the cloud.

##  IV Results based on the Gross-Pitaevskii equation

To connect more directly to our numerical GP simulations it is useful to
consider the semi-classical approximation of the quantum theory in Eq. S1
based on the GP equation:

$\displaystyle\left[J\cos(i{\partial_{y}d})-{\frac{\hbar^{2}\partial^{2}_{x}}{2m}}+g|\psi(\mathbf{r},t)|^{2}\right]\psi(\mathbf{r},t)=i{\hbar}\frac{\partial}{\partial
t}\psi(\mathbf{r},t).$ (S26)

The evolutionary dynamics found in our numerical GP simulations is seen to be
consistent with the analytics presented in the previous sections. Also
consistent is the close tie to observations from cosmological models. We can
demonstrate this through a linearization of the GP equation appropriate to the
preheating stage.

In this early stage dynamics (associated with the preheating phase in
cosmological models), the physics is governed by a parametric resonance. The
initial quenched condensate with a finite kinetic energy can be modeled by the
uniform and intrinsically oscillating scalar field,
$\psi_{0}(\mathbf{r},t)=\sqrt{n_{0}}e^{-i\omega t}.$ Here $\omega$ is the
intrinsic oscillating frequency and its value is determined by solving the
corresponding GP equation, as we show below. We assume a wavefunction of the
form, $\psi(\mathbf{r},t)=\psi_{0}({\bf{r}},t)+\delta\psi(\mathbf{r},t)$. To
linearize the theory, we treat the self-interaction term perturbatively,

$\displaystyle|\psi(\mathbf{r},t)|^{2}=\left|\sqrt{\frac{N}{V}}e^{-i\omega
t}+\delta\psi\right|^{2}\simeq n_{0}+\sqrt{\frac{N}{V}}e^{-i\omega
t}\delta\psi^{*}+h.c.+O(|\delta\psi|^{2})$ (S27)

Now we write down the equation for $\delta\psi$. Using the expansion variable
$\delta\psi$, it follows that

$\left[J\cos(i{\partial_{y}d})-{\frac{\hbar^{2}\partial^{2}_{x}}{2m}}+gn_{0}+g\sqrt{\frac{N}{V}}e^{-i\omega
t}\delta\psi^{*}+g\sqrt{\frac{N}{V}}e^{i\omega
t}\delta\psi\right]\left(\sqrt{\frac{N}{V}}e^{-i\omega
t}+\delta\psi\right)=i{\hbar}\partial_{t}\delta\psi+{\hbar}\omega\sqrt{\frac{{N}}{V}}e^{-i\omega
t}.\\\ $

Since $\psi_{0}$ satisfies the GP equation, it follows that

${\hbar}\omega=J+gn_{0},$

indicating that $\omega$ is intrinsically determined by the model parameters.
Keeping only the linear terms in $\delta\psi,\delta\psi^{*}$, one finds

$\left[J\cos(i{\partial_{y}d})-{\frac{\hbar^{2}\partial^{2}_{x}}{2m}}+2gn_{0}\right]\delta\psi+gn_{0}e^{-2i\omega
t}\delta\psi^{*}=i{\hbar}\partial_{t}\delta\psi.$ (S28)

One may write the equation in the simple form,

$E_{K}(-i\bm{\nabla})\delta\psi+g(t)\delta\psi^{*}=i{\hbar}\partial_{t}\delta\psi.$
(S29)

where
$E_{K}(-i\bm{\nabla})=J\cos(i{\partial_{y}d})-{\frac{\hbar^{2}\partial^{2}_{x}}{2m}}+2gn_{0}$
is the single-particle energy with the Hartree-Fock shift and
$g(t)=gn_{0}\exp(-2i\omega t)$ is the effective time-dependent interaction.
The time-dependence originates from the fact that the condensate is
intrinsically oscillating. Later we will see how the oscillating $g(t)$ leads
to the parametric instability. Note that in Ref. Chatrchyan  _et al._ , 2021
the modulation of the coupling (to be time-dependent) is imposed externally
and the oscillating frequency serves as a free parameter of their model. As a
contrast, the frequency $\omega$ is intrinsic in our model and has a clear
physical origin, the oscillating frequency of the condensate.

Expanding
$\delta\psi=\delta\psi(\mathbf{r},t)=\sum_{\mathbf{k},\nu}e^{i\mathbf{k}\cdot\mathbf{r}-i\nu
t}\psi_{\mathbf{k}}(\nu)$, Eq. IV then becomes

$\left[J\cos({k_{y}d})+{\frac{\hbar^{2}k_{x}^{2}}{2m}}+2gn_{0}\right]\psi_{\mathbf{k}}(\nu)+gn_{0}\psi^{*}_{-\mathbf{k}}(-\nu+2\omega)={\hbar}\nu\psi_{\mathbf{k}}(\nu).$
(S30)

With a change of variable $\nu=\nu^{\prime}+\omega$ we arrive at a matrix
equation of the form

$\left\\{\left[J\cos({k_{y}d})+{\frac{\hbar^{2}k_{x}^{2}}{2m}}+2gn_{0}-\omega\right]\hat{\sigma}_{z}+gn_{0}\hat{\sigma}_{z}\hat{\sigma}_{x}\right\\}\begin{bmatrix}\psi_{\mathbf{k}}(\nu^{\prime}+\omega)\\\
\psi^{*}_{-\mathbf{k}}(\omega-\nu^{\prime})\end{bmatrix}=\nu^{\prime}\begin{bmatrix}\psi_{\mathbf{k}}(\nu^{\prime}+\omega)\\\
\psi^{*}_{-\mathbf{k}}(\omega-\nu^{\prime})\end{bmatrix}.$ (S31)

Here $\hat{\sigma}_{x,y,z}$ is the Pauli matrix. Solving this equation, one
finds

$\nu{{}^{\prime}}(\mathbf{k})={\pm}\sqrt{\left[J\cos({k_{y}d})+{\frac{\hbar^{2}k_{x}^{2}}{2m}}+2gn_{0}-\omega\right]^{2}-(gn_{0})^{2}}.$
(S32)

When $\nu^{\prime}$ assumes imaginary values, this corresponds to the
parametric resonance condition in cosmological models. When $\omega$ is taken
to be $(J+gn_{0}){/\hbar}$, this is consistent with Eq. S6.

In this way, we have shown the existence of a parametric instability in an
isolated quantum system when the condensate has a finite intrinsic oscillating
frequency. The physical process behind the parametric resonance is that the
energy of the condensate transfers to the non-condensate particles in an
exponentially amplified way.

###  IV.1 GP derivation of the cloud thermalization stage: Comparison with the
quantum Boltzmann equation

It is instructive to use the machinery of the GP equation to characterize how
equilibration occurs. We will begin with a stage in which there is no ‘global’
phase coherence and write
$\psi(\mathbf{r},t)=\sum_{\bf{k}}\psi_{\mathbf{k}}(t)e^{i\mathbf{k}\cdot\mathbf{r}}$,
the corresponding equation of motion in terms of the complex scalar field
$\psi_{\mathbf{k}}(t)$ is then

$i\hbar\frac{\partial\psi_{\mathbf{k}}(t)}{\partial
t}=\epsilon_{\mathbf{k}}+g\sum_{\mathbf{k_{1}},\mathbf{k_{3}}}\psi_{\mathbf{k_{1}}}^{*}\psi_{\mathbf{k_{3}}}\psi_{\mathbf{k_{4}}}$
(S33)

where $\mathbf{k}_{4}=\mathbf{k}+\mathbf{k_{1}}-\mathbf{k_{3}}.$ For
simplicity, we put
$\psi_{\mathbf{k}}(t)=\tilde{\psi}_{\mathbf{k}}(t)e^{-i\epsilon_{\mathbf{k}}t/\hbar}$.
We then arrive at

$\displaystyle i\hbar\frac{\partial\tilde{\psi}_{\mathbf{k}}(t)}{\partial
t}=g$ $\displaystyle\sum
e^{i(\epsilon_{i}-\epsilon_{f})t/\hbar}\tilde{\psi}_{\mathbf{k_{1}}}^{*}\tilde{\psi}_{\mathbf{k_{3}}}\tilde{\psi}_{\mathbf{k_{4}}}$
(S34)
$\displaystyle\frac{\partial\left(\tilde{\psi}_{\mathbf{k}}^{*}\tilde{\psi}_{\mathbf{k}}\right)}{\partial
t}=$ $\displaystyle\frac{-i}{\hbar}\tilde{\psi}_{\mathbf{k}}^{*}g\sum
e^{i(\epsilon_{i}-\epsilon_{f})t/\hbar}\tilde{\psi}_{\mathbf{k_{1}}}^{*}\tilde{\psi}_{\mathbf{k_{3}}}\tilde{\psi}_{\mathbf{k_{4}}}+h.c.$

where
$\epsilon_{i}=\epsilon_{\mathbf{k_{3}}}+\epsilon_{\mathbf{k_{4}}},\,\epsilon_{f}=\epsilon_{\mathbf{k}}+\epsilon_{\mathbf{k_{1}}},\,n_{\mathbf{k}}=\psi_{\mathbf{k}}^{*}\psi_{\mathbf{k}}=\tilde{\psi}_{\mathbf{k}}^{*}\tilde{\psi}_{\mathbf{k}}.$
We now assume that in addition to the absence of ‘global’ phase coherence, the
interaction energy is quite small. As a result, the time change of
$\tilde{\psi}_{\mathbf{k}}(t)$ can be approximated to the first order in $g$,
and the change rate is much smaller than
$\epsilon_{i}/\hbar,\epsilon_{f}/\hbar$. We then obtain

$\displaystyle\tilde{\psi}_{\mathbf{k}}(t)\approx$
$\displaystyle\tilde{\psi}_{\mathbf{k}}(0)-i\pi
g\sum\delta(\epsilon_{i}-\epsilon_{f})\tilde{\psi}_{\mathbf{k_{1}}}^{*}\tilde{\psi}_{\mathbf{k_{3}}}\tilde{\psi}_{\mathbf{k_{4}}}$
(S35) $\displaystyle\frac{\partial n_{\mathbf{k}}}{\partial t}=$
$\displaystyle\frac{-ig}{\hbar}\left[\tilde{\psi}_{\mathbf{k}}^{*}(0)+i\pi
g\sum^{\prime}\delta(\epsilon_{i}^{\prime}-\epsilon_{f}^{\prime})\tilde{\psi}_{\mathbf{k_{1}^{\prime}}}\tilde{\psi}_{\mathbf{k_{3}^{\prime}}}^{*}\tilde{\psi}_{\mathbf{k_{4}^{\prime}}}^{*}\right]\sum
e^{i(\epsilon_{i}-\epsilon_{f})t/\hbar}\times\left[\tilde{\psi}_{\mathbf{k_{1}}}^{*}(0)+i\pi
g\sum^{\prime}\delta(\epsilon_{i}^{\prime}-\epsilon_{f}^{\prime})\tilde{\psi}_{\mathbf{k^{\prime}}}\tilde{\psi}_{\mathbf{k_{3}^{\prime}}}^{*}\tilde{\psi}_{\mathbf{k_{4}^{\prime}}}^{*}\right]$
$\displaystyle\times\left[\tilde{\psi}_{\mathbf{k_{3}}}(0)-i\pi
g\sum^{\prime}\delta(\epsilon_{i}^{\prime}-\epsilon_{f}^{\prime})\tilde{\psi}_{\mathbf{k_{4}^{\prime}}}^{*}\tilde{\psi}_{\mathbf{k^{\prime}}}\tilde{\psi}_{\mathbf{k_{1}^{\prime}}}\right]\times\left[\tilde{\psi}_{\mathbf{k_{4}}}(0)-i\pi
g\sum^{\prime}\delta(\epsilon_{i}^{\prime}-\epsilon_{f}^{\prime})\tilde{\psi}_{\mathbf{k_{3}^{\prime}}}^{*}\tilde{\psi}_{\mathbf{k^{\prime}}}\tilde{\psi}_{\mathbf{k_{1}^{\prime}}}\right]+h.c.$

where again both momentum and energy conservation are maintained. Now because
of the lack of phase coherence, different $\tilde{\psi}$ amplitudes have
uncorrelated random phases. Only terms with both $\psi$ and $\psi^{*}$ of the
same momentum will survive in the equation. Therefore, we find

$\displaystyle\frac{\partial n_{\mathbf{k}}}{\partial t}=$
$\displaystyle\frac{-ig}{\hbar}\left[\sum
2|\tilde{\psi}_{\mathbf{k}}|^{2}|\tilde{\psi}_{\mathbf{k_{1}}}|^{2}\right.+2\pi
ig\sum|\tilde{\psi}_{\mathbf{k_{1}}}|^{2}|\tilde{\psi}_{\mathbf{k_{3}}}|^{2}|\tilde{\psi}_{\mathbf{k_{4}}}|^{2}\delta(\epsilon_{i}-\epsilon_{f})+2\pi
ig\sum|\tilde{\psi}_{\mathbf{k}}|^{2}|\tilde{\psi}_{\mathbf{k_{3}}}|^{2}|\tilde{\psi}_{\mathbf{k_{4}}}|^{2}\delta(\epsilon_{i}-\epsilon_{f})$
(S36) $\displaystyle-2\pi
ig\sum|\tilde{\psi}_{\mathbf{k}}|^{2}|\tilde{\psi}_{\mathbf{k_{1}}}|^{2}|\tilde{\psi}_{\mathbf{k_{4}}}|^{2}\delta(\epsilon_{i}-\epsilon_{f})\left.-2\pi
ig\sum|\tilde{\psi}_{\mathbf{k}}|^{2}|\tilde{\psi}_{\mathbf{k_{1}}}|^{2}|\tilde{\psi}_{\mathbf{k_{3}}}|^{2}\delta(\epsilon_{i}-\epsilon_{f})\right]+h.c.$
$\displaystyle=\frac{4\pi
g^{2}}{\hbar}\left(n_{\mathbf{k_{1}}}n_{\mathbf{k_{3}}}n_{\mathbf{k_{4}}}+n_{\mathbf{k}}n_{\mathbf{k_{3}}}n_{\mathbf{k_{4}}}-n_{\mathbf{k}}n_{\mathbf{k_{1}}}n_{\mathbf{k_{3}}}-n_{\mathbf{k}}n_{\mathbf{k_{1}}}n_{\mathbf{k_{4}}}\right),$

where higher orders in $g$ are discarded. Comparing this result with Eq.
(S20), we find that the term
$n_{\mathbf{k_{3}}}n_{\mathbf{k_{4}}}-n_{\mathbf{k}}n_{\mathbf{k_{1}}}$ is
missing in this approximate GP approach. Thus, the GP approach works best in
the early stages of thermalization where the leading contribution is $\propto
n^{3}$. We have chosen a large particle density for our simulations to ensure
this dominance of $n^{3}$ in the early stage. We infer that there may be some
discrepancies at very late times.

##  V numerical parameters in the Gross-Pitaevskii simulations

The parameters we used for the simulations are: $V_{1}=12E_{R},V_{2}=2E_{R}$,
$gn_{0}=0.015E_{R}$, $h/T=4.75E_{R}$ where $E_{R}=\hbar^{2}(\pi/d)^{2}/2m$ is
the recoil energy unit. $d$ is the period of the lattice, which is discretized
into 64 grids in the $y$-direction with a total length of $256d$. For the
$x$-direction which is free and thus insensitive to grid resolution, we use
$64d$ for the total length to speed up the simulations. Periodic boundary
conditions are imposed for all directions. We run the dynamics of the system
for a sufficiently long time until we find the results are stable. The longest
running time is 6000T.

Additionally, random noise terms have been added to the system to simulate
quantum fluctuations. Our simulations involve Graphics Processing Unit-based
quasispectral, split-step method to solve the GP equation based on fast
Fourier transforms. More details can be found in “Pseudo-spectral solution of
nonlinear Schroedinger equations” in the Journal of Computational Physics 87,
pages 108-125 (1990).

## References

  * Feng  _et al._ (2018) L. Feng, L. W. Clark, A. Gaj, and C. Chin, Nature Physics 14, 269 (2018).
  * Note (1) To observe why $\mathbf{q}_{1}=-\mathbf{q}_{2}$ is the leading contribution, one may define $f(Q)$ as the multiplicity for the space $q_{1}+q_{2}=Q$, i.e., $f(Q)=\text{Number of elements in }\\{(q_{1},q_{2})|\,q_{1}+q_{2}=Q,\,q_{1,2}\in\Lambda\\}.$ One can check that, for the ring-shape $\Lambda$, $f(Q)$ is on the order of $1$ if $Q\neq 0$ while for $Q=0$, $f(Q)$ is of the order $L\gg 1$ .
  * Chatrchyan  _et al._ (2021) A. Chatrchyan, K. T. Geier, M. K. Oberthaler, J. Berges, and P. Hauke, Phys. Rev. A 104, 023302 (2021).