# Opinion Dynamics in Financial Markets via Random Networks

Mateus F. B. Granha<EMAIL_ADDRESS>Física de Materiais, Universidade de
Pernambuco, Recife, PE 50720-001, Brazil  André L. M. Vilela
<EMAIL_ADDRESS>Universidade de Pernambuco, Recife, PE 50720-001, Brazil
Center for Polymer Studies and Department of Physics, Boston University,
Boston, MA 02215, USA  Chao Wang<EMAIL_ADDRESS>College of Economics
and Management, Beijing University of Technology, Beijing, 100124, China
Kenric P. Nelson<EMAIL_ADDRESS>Photrek LLC, Watertown, MA 02472, USA  H.
Eugene Stanley<EMAIL_ADDRESS>Center for Polymer Studies and Department of
Physics, Boston University, Boston, MA 02215, USA

###### Abstract

We investigate the financial market dynamics by introducing a heterogeneous
agent-based opinion formation model. In this work, we organize the individuals
in a financial market by their trading strategy, namely noise traders and
fundamentalists. The opinion of a local majority compels the market exchanging
behavior of noise traders, whereas the global behavior of the market
influences the fundamentalist agents’ decisions. We introduce a noise
parameter $q$ to represent a level of anxiety and perceived uncertainty
regarding the market behavior, enabling the possibility for an adrift
financial action. We place the individuals as nodes in an Erdös-Rényi random
graph, where the links represent their social interaction. At a given time,
they assume one of two possible opinion states $\pm 1$ regarding buying or
selling an asset. The model exhibits such fundamental qualitative and
quantitative real-world market features as the distribution of logarithmic
returns with fat-tails, clustered volatility, and long-term correlation of
returns. We use Student’s t distributions to fit the histograms of logarithmic
returns, showing the gradual shift from a leptokurtic to a mesokurtic regime,
depending on the fraction of fundamentalist agents. We also compare our
results with the distribution of logarithmic returns of several real-world
financial indices.

Econophysics, Sociophysics, Monte Carlo simulation, Phase transitions, Complex
Networks

###### pacs:

89.65.Gh Economics; econophysics, financial markets, business and management,
87.23.Ge Dynamics of social systems, 05.10.Ln Monte Carlo methods, 64.60.Cn
Order-disorder transformations

††preprint: APS/123-QED

##  I INTRODUCTION

For decades now, especially over the last years, financial markets and
economic systems have fascinated researchers and investors worldwide. The
effectiveness of methods and techniques from Statistical Mechanics has been a
compelling element for the comprehension of financial dynamics as a complex
system stanley2000introduction ; Stanley2002 ; Farmer2005 ; Baldovin2007 ;
Li2020 . Simultaneously, it has been responsible for developing the
interdisciplinary field of Econophysics, where agent-based models have a
significant contribution. In this framework, several models use a set of
elementary rules to represent the agent behavior and interactions. Thus,
yielding intense emergent collective phenomena bonabeau2002agent ;
macal2005tutorial ; macal2009agent ; vilela2019majority ; zubillaga2019three ;
Feng2012 ; Zhao2011 .

Oliveira de1992isotropic ; Oliveira1993  pioneered in the study of opinion
dynamics modeled by the majority-vote model (MVM) with noise, a sociological
adaptation of a magnetic spin model that presents similar critical behavior.
Lima et al., Campos et al. and Pereira et al. lima2008majority ;
campos2003small ; pereira2005majority  investigated the effects of complex
networks on the MVM social dynamics, thus revealing an ordered phase over the
increasing complexity of these networks. Vilela et al. vilela2019majority
incorporated the rational and emotional behavior of agents in financial
markets introducing the global-vote model, inspired by the majority-vote
dynamics.

The MVM consists of a comprehensive approach to the study of social dynamics
with interactions embedded in regular and complex networks de1992isotropic ;
Oliveira1993 ; santos1995anisotropic ; vieira2016phase ; pereira2005majority ;
lima2008majority ; campos2003small ; vilela2009 ; vilela2018effect ;
vilela2021 . In this agent-based model with two states, the individual’s
opinion in a given time instant may assume one of the two values, $\pm 1$,
regarding some social discussion. Furthermore, an agent in the social network
assumes the opinion of the majority of its neighboring spins with probability
$(1-q)$ and the opposite opinion with probability $q$. The variable $q$ stands
for the noise parameter of the model and measures the social unrest or social
temperature of the system. The MVM exhibits a second-order phase transition
for a critical noise value $q=q_{c}\approx 0.075$ in square lattice networks
of social interactions de1992isotropic .

Complex Networks consist of a natural preference to the study of real-world
complex systems, such as climate analysis, biological neural connections,
World Wide Web, public transportation, airline networks, financial markets,
among others feldhoff2015complex ; reijneveld2007application ;
barabasi2000scale ; an2014synchronization . In this context, the random graph
networks describe the topology of a set of $N$ nodes connected with chance
$p$, thus introducing a fundamental probability distribution over graphs. The
Erdös-Rényi algorithm is one well-known method for the assembly of random
graph networks. The method connects an initial set of $N$ isolated nodes by
adding a total of $pN(N-1)/2$ random links between them, while double
connections are forbidden erdos1960evolution ; bollobas1985random .

In Figure 1, we illustrate the Erdös-Rényi process of building a random graph
network with $N=10$ and $\left<k\right>=3$. We also present the averaged
degree distribution $P(k)$ over ten networks with $N=2\times 10^{4}$ nodes and
several values of $\left<k\right>$. The lines correspond to Poisson fits for
the data with average degree $\left<k\right>=pN$, in agreement for the degree
distribution of random graphs for large networks.

Figure 1: Illustration of the Erdös-Rényi process of building random graph
networks. From (a) to (d), $N=10$ isolated nodes are connected by adding a
total of $pN(N-1)/2$ random links between them. While double connections are
forbidden, the final network has $\left<k\right>=3$. In (e) we show the degree
distribution $P(k)$ averaged over ten networks for $\left<k\right>=6,8,10$ and
$20$ with $N=2\times 10^{4}$. The lines represent Poisson fits for the data.

Modern financial markets operate as a complex system. Its dynamics depend not
only on the rational strategy of the individual but also on the emotional
circumstances that define the investor’s psyche. At first, individual
decisions might seem challenging to model; still, social agents tend to follow
a herding behavior, as they feel sheltered when the crowd is endorsing their
choices. In this way, trader social dynamics is reasonably attainable
cont2000herd ; sznajd2002simple . This bias to follow a significant group is
often embraced by less experienced agents, denominated noise traders. They are
highly susceptible to the dominant opinion and tend to overreact to fresh news
regarding buying or selling.

In contrast, other traders, denominated noise contrarian traders or
fundamentalists, follow the global minority in a given market as an investment
strategy. Thus they buy when stock prices decline and sell when the prices
increase sznajd2002simple ; lux1999scaling ; bornholdt2001expectation ;
kaizoji2002dynamics ; takaishi2005simulations . Therefore, their decision-
making process drives the asset price to its fundamental value.

Inspired by the global-vote model for financial markets vilela2019majority ,
we propose to investigate the influence of random graph networks on the time
evolution of financial markets quantity measures. This work employs an opinion
formation model to study the economic mechanics of noise traders and
contrarian individuals interacting through a random graph of social influence.

This work is organized as follows. In Section II, we describe the global-vote
opinion formation model for financial markets and present the relevant
quantities analyzed in our simulations. In Section III, we present the
numerical results obtained along with the corresponding discussions. Section
IV concludes with our final remarks.

##  II THE MODEL

We represent the financial agents on the market and their interactions using a
complex network structure. We place $N$ agents on the nodes of a random graph
network, where the links represent the interaction between the agents in the
market. We map the agent’s financial decision (opinion or option) at a given
time $t$ by a spin variable, which may assume one of two values: $+1$ or $-1$,
regarding buying or selling an asset. Furthermore, to model the essential
dynamics of real-world financial markets, we randomly distribute two sets of
individuals: a fraction $f$ of fundamentalist traders, also referred to as
noise contrarians, and the remaining fraction $1-f$ of noise traders. We use
the spin variables $\lambda$ and $\alpha$ to stand for the financial option of
noise traders and contrarian agents, respectively.

We introduce a noise parameter variable $q$ to model a financial level of
anxiety and perceived uncertainty present in the market. In this study, $q$
influences contrarian and noise traders’ decisions and represents the
probability of performing an inaccurate financial action. In other words, $q$
is the chance for an agent to choose the opposite of his standard strategy
when negotiating in the financial market.

###  II.1 Noise Traders

A noise trader individual $i$, with opinion $\lambda_{i}$, assumes the same
opinion of the majority of its neighbors with probability $1-q$, and the
opposite option with probability $q$. We write the opinion flipping
probability for the noise trader agent $i$ as follows

$\omega(\lambda_{i})=\frac{1}{2}\left[1-(1-2q)c_{\lambda}\lambda_{i}\textrm{sgn}(m)\right],$
(1)

where $\textrm{sgn}(x)=-1,0,+1$ for $x<0$, $x=0$ and $x>0$, respectively.
Here, $c_{\lambda}$ stands for the agent strategy, where we set
$c_{\lambda}=+1$ to model the noise traders trend to agree with the local
majority of their neighbors. The variable $m$ quantifies the local predominant
opinion, or local magnetization, defined as

$m=\sum_{\delta=1}^{k_{i}}\lambda_{i+\delta}.$ (2)

The summation is over all the $k_{i}$ neighboring agents connected to the
trader at node $i$, with opinion $\lambda_{i}$. From Eq. (1), when $q=0$ the
noise trader adopts its local predominant opinion. When we increase the market
anxiety parameter $q$, the noise trader tends to follow the opposite opinion
of its local majority.

###  II.2 Noise Contrarians

A fundamentalist agent tends to follow the market’s minority opinion with
probability $1-q$, while following the majority opinion with chance $q$. We
define the prevailing option of the system as a global magnetization, which
influences noise contrarian traders’ financial option. The option of a noise
contrarian agent $j$, $\alpha_{j}$, flips with probability

$\omega(\alpha_{j})=\frac{1}{2}\left[1-(1-2q)c_{\alpha}\alpha_{j}\textrm{sgn}(M)\right],$
(3)

where $c_{\alpha}$ is defined as the fundamentalist’s strategy, and since
these individuals tend to agree with the global minority of the system,
$c_{\alpha}=-1$, for all noise contrarian agents.

The variable $M$ measures the average market opinion of the system, thus
revealing the economic order. The magnetization $M$ accounts for the financial
opinion of every agent on the market, and is evaluated as

$M=\frac{1}{(N_{\lambda}+N_{\alpha})}\left(\sum_{i=1}^{N_{\lambda}}\lambda_{i}+\sum_{j=1}^{N_{\alpha}}\alpha_{j}\right),$
(4)

where $N_{\lambda}=N(1-f)$ stands for the number of noise trader agents in the
network and $N_{\alpha}=Nf$ represents the number of noise contrarian agents.
For $M=1$, every agent on the market has an option equal to $+1$. Similarly,
$M=-1$ denotes a market configuration where every agent opinion is equal to
$-1$. $M=0$ represents the case where half of the agents have opinion $+1$ and
the other half $-1$. Additionally, intermediate values of $M$ infer the
dominant opinion of the market.

We remark that when $q=0$, a contrarian agent always follow the opposite
opinion of the global magnetization: buy (sell) when the majority sells (buy).
As $q$ increases, there is a greater probability for the contrarian agent to
follow the majority, adopting the opposite of his inherent financial market
strategy.

In the global-vote dynamics for financial markets, we recover the standard
majority-vote model on random graphs when the fraction of contrarians $f$ is
zero. In this case, we observe a order-disorder phase transition at a critical
point $q_{c}$ for each value of the average connectivity $\left<k\right>$ in
the thermodynamic limit $N\rightarrow\infty$ pereira2005majority . This $f=0$
system presents an ordered phase, with large clusters of agents that share the
same opinion for values of noise $q$ below the critical point $q_{c}$. By
increasing $q$, the same opinion clusters fade, and the magnetization (order
parameter) approaches zero.

In this work, we analyze the global-vote dynamics for financial markets on
random graphs for several values of the average connectivity $\left<k\right>$
and noise contrarian traders fraction $f\neq 0$. We set $q$ near
$q_{c}(f=0,\left<k\right>)$ to model real-world market dynamics adequately
since previous investigations suggest the strong market phase emerge when the
system is close to its critical melting point bornholdt2001expectation ;
vilela2019majority .

##  III RESULTS AND DISCUSSION

We perform Monte Carlo simulations on Erdös-Rényi random graphs of size
$N=10201$ and several values of the average connectivity $\left<k\right>$. For
the network of financial interactions, we randomly place a fraction $1-f$ of
noise trader on the network nodes, and we occupy the remaining fraction $f$
with noise contrarian agents. In this way, $N=N_{\lambda}+N_{\alpha}$.

This work considers that a real-world market presents a small fraction of
noise contrarian agents. To this extent, we investigate the influence of the
noise contrarian agents for $q$ near its critical value for the $f=0$ case,
obtained in previous studies within several values of $\left<k\right>$
pereira2005majority , where our model exhibits key real-world market features
vilela2019majority . Additionally, we consider small values for the average
number of connections of an agent in a market, i.e., small values of
$\left<k\right>$, supported by the agreement between the data of real-world
markets and the results of our simulations.

The market dynamics evolve as follows. We select a randomly chosen agent and
update its opinion accordingly with the probabilities given by Eqs. (1) or (3)
for a noise trader or a noise contrarian agent, respectively. We repeat this
process $N$ times, thus defining a unity of time in one Monte Carlo step
(MCS). This way, each agent’s opinion is updated once on average on one MCS.
To discard the transient regime, we allow this dynamic to run during $10^{3}$
MCS, and we perform our analysis in the subsequent $10^{5}$ MCS.

In the financial context, we shall relate the global magnetization of the
system to the aggregate excess demand of a particular asset, for its impact on
stock prices bornholdt2001expectation ; kaizoji2002dynamics ;
takaishi2005simulations . Thus a positive demand ($M>0$) causes prices to
rise, while a negative demand ($M<0$) causes prices to fall. Moreover, markets
that fluctuate around equilibrium cont2000herd  exhibit an average excess
demand that oscillates around zero.

The financial return time series represents the price variation of a given
asset over time. Positive (negative) returns relate to profit (loss) during
the period of analysis. We quantify the logarithmic return at a given time $t$
as follows

$r(t)=\textrm{log}\left[\left|M(t)\right|\right]-\textrm{log}\left[\left|M(t-1)\right|\right].$
(5)

In Figure 2, we present the logarithmic return for the closing values of the
S&P 500 index in US dollars, from Nov 01, 2012 to Dec 08, 2020. Similar to any
other traditional asset, the S&P 500 index depends mainly on market supply and
demand, a fundamental financial mechanism that yields several periods of
strong return variations. Periods of significant fluctuations of returns are
compressed in time, denoting the clustered volatility effect for the analyzed
index price. This financial phenomenon is known as volatility clustering, and
it can be elucidated by Mandelbrot’s observation that “large changes tend to
be followed by large changes - of either sign - and small changes tend to be
followed by small changes” mandelbrot1963variation .

Figure 2 also shows that the period around t = 1850 presents the most
considerable time-series return volatility. It denotes the Coronavirus
pandemic phase, where government officials halted the economic activity. The
panic and uncertainty triggered by the financial impacts of such measures led
to a expressive stock market crash vega2021 ; wei2021 .

Figure 2: Logarithmic returns for the daily closing price of S&P 500 index in
US dollars from Nov 01, 2012 to Dec 08, 2020. The high volatility observed for
the period around t = 1850 represents the Coronavirus Stock Market Crash of
2020.

In Figure 3, we analyze the influence of the average connectivity
$\left<k\right>$ and the noise contrarian fraction $f$ on the market order of
the model on random graphs. Figure 3(a) exhibits two distinct market phases: a
“strong market phase” for $f=0.20$ (dark blue), where the system is highly
volatile and the magnetization exhibits a irregular wave pattern; and a “weak
market phase” for $f=0.70$ (yellow), where magnetization values are roughly
randomly distributed vilela2019majority . This result demonstrates that the
increase of contrarians tends to stabilize the market dynamics, where the high
volatile periods are present but to a substantially lower extent.
Additionally, the increase of the agent’s average number of connections
$\left<k\right>$ also yields stochastic patterns for the market demand, but
driving a contraction on the amplitude of the fluctuations, as observed in
Figure 3(b).

Figure 3: Time series of the order parameter $M$ for two sets of parameters:
(a) $\left<k\right>=6$ for $q=0.240$, and (b) $\left<k\right>=50$ for
$q=0.411$.

Figure 4, displays the logarithmic returns of the absolute value of the order
parameter. In financial markets, we define volatility as a measure of the
variation around the average returns observed in an asset’s time series. In
this figure, we note intensive market fluctuations for $f=0.20$, embodied by
the large spikes in the plot. The result also shows clustered volatility, a
real-world market feature. We observe that the presence of contrarians tends
to stabilize the market, indicated by attenuated fluctuations of the returns.
Furthermore, when comparing Figures 4(a) and 4(b), it becomes clear that
increasing the average connectivity of the network will similarly increase the
market’s volatility and simultaneously deviate the behavior from the expected
of real-world financial markets. Figure 4 also shows that periods of high
volatility tend to be clustered, for lower values of $f$.

Figure 4: Logarithmic returns of the absolute magnetization for (a)
$\left<k\right>=6$ and $q=0.240$, (b) $\left<k\right>=50$ and $q=0.411$.

In order to quantify the effects of volatility clustering, we shall compute
the autocorrelation of absolute returns as follows

$A(\tau)=\frac{\sum_{t=\tau+1}^{T}\left[\left|r(t)\right|-|\overset{-}{r}|\right]\left[\left|r(t-\tau)\right|-|\overset{-}{r}|\right]}{\sum_{t=1}^{T}\left[\left|r(t)\right|-|\overset{-}{r}|\right]^{2}},$
(6)

where $1\leq\tau\leq 10^{5}\ \textrm{MCS}$ is the time-step difference between
observations, $T=10^{5}\ \textrm{MCS}$ is the time of simulation, $r(t)$ is
the return at a time $t$ and $\overset{-}{r}$ the average value of the return.
The function defined by Eq. (6) measures nonlinear correlations in
observations of the absolute value of log-returns as a function of the time
delay between them. Figure 5 displays the autocorrelation for several values
of $f$, and it shows that the returns present a strong correlation in time
with exponential decay kaizoji2002dynamics ; takaishi2005simulations ;
bornholdt2001expectation ; vilela2019majority .

Figure 5: Linear-log plot of the autocorrelation of absolute logarithmic
returns for (a) $\left<k\right>=6$ and $q=0.240$, and (b) $\left<k\right>=8$
and $q=0.275$. The dashed red lines represent exponential fits for the data.
Figure 6: Linear-log plot of the autocorrelation of absolute logarithmic
returns of the closing values for several financial indices. Also shown is the
autocorrelation of absolute log-returns for $\left<k\right>=6$, $q=0.240$ and
$f=0.20$ (dark blue).

To illustrate this exponential behavior, we perform an exponential fit of the
autocorrelation function $A(q,f,t)\sim\textrm{exp}(-t/t_{0})$ for $f=0.20$,
and we obtain $1/t_{0}\approx 2.5\times 10^{-7}$ and $1/t_{0}\approx 3.0\times
10^{-7}$ for Figures 5(a) and 5(b), respectively. Other values for the
fraction of noise contrarians yield similar results.

In Figure 6, we compare the results of our model with real-world financial
indices. We calculate the autocorrelation function for the daily log-returns
of the closing values of the indices: Dow Jones (DJI) observed from Jan 29,
1985 to Dec 02, 2020; Ibovespa (BVSP) observed from July 24, 1993 to Dec 02,
2020; Nikkei (N225) observed from Jan 05, 1980 to Dec 02, 2020; S&P 500 (GSPC)
observed from Jan 02, 1985 to Dec 02, 2020; Nasdaq (IXIC) observed from Oct
01, 1985 to Dec 02, 2020. We analyze each index for roughly $7\times 10^{3}$
days. We note that the returns for each investigated index display an
expressive correlation which decays exponentially in time.

Figure 7 displays the histogram of the log-returns for two sets of parameters
and several values of the noise contrarian fraction $f$ in $10^{5}$ MCS. The
real-world market systems display fat-tailed distributions as a reflection of
lower non-zero probabilities of obtaining above or below-average returns. We
observe such behavior only for lower values of the fraction of contrarians
$f$, substantiating a real-world market region for our agent-based dynamics.
To quantify the return distributions, we perform a statistical analysis of the
data. We obtain that the kurtosis $\textrm{K}(\left<k\right>,f)$ for $f=0.20$
(strong market phase) are $\textrm{K}(6,0.20)=5.85$ and
$\textrm{K}(8,0.20)=4.98$. For $f=0.70$ (weak market phase),
$\textrm{K}(6,0.70)=2.41$ and $\textrm{K}(8,0.70)=2.37$. We remark that for a
normal distribution $\textrm{K}=3$, thus, increasing the fraction of
contrarians $f$ gradually shifts the system’s behavior, from a leptokurtic
(fat-tailed) regime, to a mesokurtic (Gaussian) regime.

To qualify the distributions of Figure 7, we perform a comparative normal
quantile-quantile (Q-Q) plot. Fig. 8 displays Q-Q plots for the distribution
of log-returns using several values for the fraction of fundamentalists $f$.
The red line represents the expected results of a Gaussian distribution, and
if a particular distribution exhibits similar behavior, its data points should
lie on that reference line. For lower values of $f$, for instance $f\leq
0.40$, both Figures 8(a) and 8(b) display a non-linear behavior in the Q-Q
plot, thus indicating that the distributions of log-returns feature fat-tails,
and therefore, non-Gaussian. In such distributions, one observes a greater
chance of obtaining high return values than the expected for normal
distributions, consistent with real-world market systems
stanley2000introduction .

Figure 7: Distribution of logarithmic returns for $10^{5}$ MCS and several
values of the fraction of contrarians $f$: (a) $\left<k\right>=6$ and
$q=0.240$, and (b) $\left<k\right>=8$ and $q=0.275$. Figure 8: Normal
quantile-quantile plots of the logarithmic return distributions for (a)
$\left<k\right>=6$ and $q=0.240$, and (b) $\left<k\right>=8$ and $q=0.275$.
Here we use $10^{5}$ MCS. Figure 9: Plot of the absolute log-return
distributions $F(q,\ f,\ r)$ for several values of the fraction $f$ and
$10^{5}$ MCS . The lines for $f\leq 0.40$ correspond to Student’s t fits,
whereas the lines for $f\geq 0.5$ correspond to Gaussian fits. We use (a)
$\left<k\right>=6$ and $q=0.240$, and (b) $\left<k\right>=8$ with $q=0.275$.

Figure 9 displays the probability distribution of the absolute log-returns for
several values of $f$ in $10^{5}$ MCS with $\left<k\right>=6$ and
$\left<k\right>=8$. We find that Student’s t-distributions can depict the
strong market’s phase data ($f\leq 0.40$), while the weak market’s phase is
well fitted by Gaussian distributions vilela2019majority . Assuming that the
mean values $\mu$ of the distributions are zero, as illustrated by Figure 10,
we perform generalized Student’s t fits for the data
$S(r;\nu,\mu,\sigma)=S(r;\nu,\sigma)$, where $\nu$ is the degree of freedom
and $\sigma$ represents the scale. The function $S$ is defined as follows

$S(r;\nu,\sigma)=\frac{1}{\sqrt{\nu\sigma^{2}}\textrm{B}\left(\frac{\nu}{2},\frac{1}{2}\right)}\left(1+\frac{r^{2}}{\nu\sigma^{2}}\right)^{-\frac{\nu+1}{2}},$
(7)

where

$\textrm{B}\left(\frac{\nu}{2},\frac{1}{2}\right)=\int_{0}^{1}t^{\frac{\nu}{2}-1}(1-t)^{-\frac{1}{2}}\textrm{dt},$
(8)

is the Beta function. Values obtained for the degree $\nu$ and scale $\sigma$
are displayed in Table 1. As the fraction of contrarians increase, we find
that the absolute log-return distributions are well fitted by a Gaussian
distribution $g(r)$

$g(r)=\frac{1}{\sigma\sqrt{2\pi}}\
\textrm{e}^{-\frac{1}{2}\left(\frac{r-\mu}{\sigma}\right)^{2}},$ (9)

where $\sigma$ is the standard deviation and $\mu$ is the mean, which we shall
consider as zero once more. In Table 2, we display the Gaussian fit
information.

|  |  $\left<k\right>=6$, $q=0.240$ |  |  $\left<k\right>=8$, $q=0.275$
---|---|---|---|---
$f$ |  | $\nu$ | $\sigma$ |  | $\nu$ | $\sigma$
0.20 |  | 5.0 $\pm$ 0.3 | 0.80 $\pm$ 0.01 |  | 3.0 $\pm$ 0.1 | 0.785 $\pm$ 0.008
0.25 |  | 7.4 $\pm$ 0.5 | 0.558 $\pm$ 0.006 |  | 6.3 $\pm$ 0.4 | 0.621 $\pm$ 0.006
0.30 |  | 12 $\pm$ 2 | 0.444 $\pm$ 0.006 |  | 10 $\pm$ 1 | 0.559 $\pm$ 0.008
0.40 |  | 13 $\pm$ 2 | 0.252 $\pm$ 0.005 |  | 38 $\pm$ 2 | 0.319 $\pm$ 0.006
Table 1: Correlation between the scale $\sigma$ and degree $\nu$ as a function of the fraction of contrarians f, the average connectivity $\left<k\right>$ and noise $q$. |  |  $\left<k\right>=6$, $q=0.240$ |  |  $\left<k\right>=8$, $q=0.275$
---|---|---|---|---
$f$ |  | $\sigma$ |  | $\sigma$
0.50 |  | 0.1916 $\pm$ 0.0008 |  | 0.2254 $\pm$ 0.0009
0.70 |  | 0.1108 $\pm$ 0.0004 |  | 0.1284 $\pm$ 0.0005
Table 2: Dependence of the standard deviation $\sigma$ on the fraction f, on
the average connectivity $\left<k\right>$ and on noise parameter $q$.

Figure 10 displays the cumulative distribution of logarithmic returns $\Phi$
for the data in $10^{5}$ MCS, from which we imply that the mean for the
distributions remains zero for $\left<k\right>=6$ with $q=0.240$, and
$\left<k\right>=8$ with $q=0.275$ for several values of $f$.

Figure 10: Plot of the cumulative distribution $\Phi$ of log-returns in
$10^{5}$ MCS for (a) $\left<k\right>=6$ and $q=0.240$, and (b)
$\left<k\right>=8$ and $q=0.275$. In the inset we analyze the details of the
plot for values of $r(t)$ near zero.

Figure 11: Plot of the absolute log-return distributions $F(r)$ of the closing
values of several financial indices. Also shown is the absolute log-return
distributions for $\left<k\right>=6$, $q=0.240$, $f=0.20$. The lines
correspond to Student’s t fits.

In Figure 11, we show the distribution of the absolute log-returns of the
closing values of several financial indices to the results found in our
investigation. We analyze each index for roughly $7\times 10^{3}$ days and use
Student’s t-distributions to fit the data, represented by the lines in the
plots. In Table 3, we show the values for the degree $\nu$ and scale $\sigma$.
We note that the degree of freedom primarily determines the shape of each
distribution, and comparison shows that our model behaves as the expected for
real-world financial markets. On the other hand, the scale parameter measures
fluctuations around the mean of the distribution, and comparison displays a
disparity in the order of magnitude obtained in the fits. Analyzing Figures 2
and 4, we observe that the daily log-returns of financial markets are one
order of magnitude smaller than the log-returns obtained in our simulations,
and thus accountable for the divergence in the results. Hence, we observe that
our model represents satisfactorily real-world financial market behavior.

Index | $\nu$ | $\sigma$
---|---|---
Dow Jones | 9.5 $\pm$ 0.2 | 0.01248 $\pm$ 8 $\times 10^{-5}$
Ibovespa | 4.2 $\pm$ 0.3 | 0.0197 $\pm$ 6 $\times 10^{-4}$
Nikkei | 15.3 $\pm$ 0.6 | 0.0202 $\pm$ 2 $\times 10^{-4}$
S&P 500 | 15.6 $\pm$ 0.5 | 0.0151 $\pm$ 1 $\times 10^{-4}$
Nasdaq | 6.5 $\pm$ 0.3 | 0.0180 $\pm$ 3 $\times 10^{-4}$
Table 3: Values for the scale $\sigma$ and degree $\nu$ for several financial
indices with Student’s t fits.

We also perform an extensive investigation of the model for several values of
the average connectivity $\left<k\right>$ and different values of the noise
parameter $q$. Figure 12 displays a multi-plot table of the histograms of
logarithmic returns for several ($q,f,\left<k\right>$) triplets. Columns
correspond to values of $q$, below criticality $q<q_{c}(\left<k\right>)$
(left), at criticality $q=q_{c}(\left<k\right>)$ (center) and above
criticality with $q>q_{c}(\left<k\right>)$ (right). In this work, we use the
values of $q=q_{c}(\left<k\right>)$ obtained in previous investigations of the
MVM on random graphs, where the fraction of contrarian agents $f$ is set to
zero pereira2005majority . The values of noise above and below criticality are
taken to be $q_{c}(1\pm 10\%)$.

We observe that distributions exhibit fat-tails for lower values of $f$, and
especially for small values of $\left<k\right>$, i.e., $\left<k\right>=6$ and
$\left<k\right>=8$, and at its correspondent critical value of $q$. On the
contrary, such behavior is lost for higher values of the average connectivity
and values of $q$ that deviate from criticality. This result illustrates and
supports our choice for $\left<k\right>$ and $q_{c}$ used in this
investigation.

Figure 12: Distributions of logarithmic returns for several of the average
degree of connectivity $\left<k\right>$ in the vicinity of
$q_{c}(\left<k\right>)$ for several values of $f$: 0.20; 0.25; 0.30; 0.50;
0.70 (dark blue, purple, violet, pink, orange, yellow, respectively). Here, we
use $10^{5}$ MCS.

We conclude that the adoption of random graph networks on the global-vote
model for financial markets has demonstrated to be effective. Our model is
able to reproduce qualitatively and quantitatively real-world market features
for lower values of the average connectivity $\left<k\right>$ and near
criticality $q(\left<k\right>)=q_{c}(\left<k\right>)$.

We remark that other combination of values for $\left<k\right>$, $f$ and $q$
might yield similar results. Nevertheless, our particular choice adopts simple
key ideas: limited and small number of interacting agents, near critical
values of the noise parameter $q$ and small number of contrarian agents in the
market. Despite the model’s simplicity, it has shown its capability of
characterizing the mechanisms that drive social behavior and decision making
in economic systems.

##  IV CONCLUSION AND FINAL REMARKS

This work proposes a generalization of the two-state global-vote model for
financial markets on random networks. The global-vote model suggests that any
stock market dynamics consist primarily of different agent strategies driven
by economic and social interactions. In its standard version, the stock market
consists of a heterogeneous population with two distinct kinds of investors:
noise traders, who follow the local majority of its neighbors, and noise
contrarian traders, influenced by the global minority of the system
vilela2019majority .

We aim to investigate the return’s distribution dependence on the average
connectivity $\left<k\right>$ of the random graph network. We relate
variations of the global magnetization of the system to the daily return of a
given asset. Our simulations reproduce the typical qualitative and
quantitative real-world financial time series. Thus, yielding key features as
fat-tailed distributions of returns, volatility clustering and long-term
memory volatility. We demonstrate that higher values for the average
connectivity $\left<k\right>$, or the noise contrarian fraction $f$, as well
as values of $q$ far from criticality, may eliminate similar real-world
behavior in our model.

This investigation combines two distinct fields of research - Econophysics and
Sociophysics - to better understand the mechanisms at play in the complexity
of the human decision-making process, which drives economic systems’ dynamics.
Consequently, it demonstrates advancements towards comprehending real-world
heterogeneous complex systems, such as financial markets.

Acknowledgements

The authors acknowledge financial support from Brazilian and Chinese
institutions and funding agents UPE, FACEPE (APQ-0565-1.05/14,
APQ-0707-1.05/14), CAPES, CNPq, National Natural Science Foundation of China
(72071006), the International Postdoctoral Exchange Fellowship Program
(20170016) and Beijing Social Science Foundation of China (16JDGLC005). The
Boston University Center for Polymer Studies is supported by NSF Grants
PHY-1505000, CMMI-1125290, and CHE-1213217, by DTRA Grant HDTRA1-14-1-0017,
and by DOE Contract DE-AC07-05Id14517.

## References

  * (1) H. E. Stanley and R. N. Mantegna, An introduction to econophysics.  Cambridge University Press, Cambridge, 2000.
  * (2) H. E. Stanley, L. A. N. Amaral, S. V. Buldyrev, P. Gopikrishnan, V. Plerou, and M. A. Salinger, “Self-organized complexity in economics and finance,” Proceedings of the National Academy of Sciences, vol. 99, no. suppl 1, pp. 2561–2565, 2002.
  * (3) J. D. Farmer, P. Patelli, and I. I. Zovko, “The predictive power of zero intelligence in financial markets,” Proceedings of the National Academy of Sciences, vol. 102, no. 6, pp. 2254–2259, 2005.
  * (4) F. Baldovin and A. L. Stella, “Scaling and efficiency determine the irreversible evolution of a market,” Proceedings of the National Academy of Sciences, vol. 104, no. 50, pp. 19741–19744, 2007.
  * (5) Y. Li, A. L. M. Vilela, and H. E. Stanley, “The institutional characteristics of multifractal spectrum of china’s stock market,” Physica A: Statistical Mechanics and its Applications, vol. 550, p. 124129, 2020.
  * (6) E. Bonabeau, “Agent-based modeling: Methods and techniques for simulating human systems,” Proceedings of the National Academy of Sciences, vol. 99, no. suppl 3, pp. 7280–7287, 2002.
  * (7) C. M. Macal and M. J. North, “Tutorial on agent-based modeling and simulation,” in Proceedings of the 2005 Winter Simulation Conference, p. 14, IEEE, 2005.
  * (8) C. M. Macal and M. J. North, “Agent-based modeling and simulation,” in Proceedings of the 2009 Winter Simulation Conference, pp. 86–98, IEEE, 2009\.
  * (9) A. L. M. Vilela, C. Wang, K. P. Nelson, and H. E. Stanley, “Majority-vote model for financial markets,” Physica A: Statistical Mechanics and its Applications, vol. 515, pp. 762–770, 2019.
  * (10) B. J. Zubillaga, A. L. M. Vilela, C. Wang, K. P. Nelson, and H. E. Stanley, “A three-state opinion formation model for financial markets,” Physica A: Statistical Mechanics and its Applications, vol. 588, p. 126527, 2022.
  * (11) L. Feng, B. Li, B. Podobnik, T. Preis, and H. E. Stanley, “Linking agent-based models and stochastic models of financial markets,” Proceedings of the National Academy of Sciences, vol. 109, no. 22, pp. 8388–8393, 2012.
  * (12) L. Zhao, G. Yang, W. Wang, Y. Chen, J. P. Huang, H. Ohashi, and H. E. Stanley, “Herd behavior in a complex adaptive system,” Proceedings of the National Academy of Sciences, vol. 108, no. 37, pp. 15058–15063, 2011.
  * (13) M. J. de Oliveira, “Isotropic majority-vote model on a square lattice,” Journal of Statistical Physics, vol. 66, no. 1-2, pp. 273–281, 1992.
  * (14) M. J. de Oliveira, J. F. F. Mendes, and M. A. Santos, “Nonequilibrium spin models with ising universal behaviour,” Journal of Physics A: Mathematical and General, vol. 26, no. 10, p. 2317, 1993.
  * (15) F. W. S. Lima, A. O. Sousa, and M. A. Sumuor, “Majority-vote on directed erdős–rényi random graphs,” Physica A: Statistical Mechanics and its Applications, vol. 387, no. 14, pp. 3503–3510, 2008.
  * (16) P. R. A. Campos, V. M. de Oliveira, and F. G. B. Moreira, “Small-world effects in the majority-vote model,” Physical Review E, vol. 67, no. 2, p. 026104, 2003.
  * (17) L. F. C. Pereira and F. G. B. Moreira, “Majority-vote model on random graphs,” Physical Review E, vol. 71, no. 1, p. 016123, 2005.
  * (18) M. A. Santos and S. Teixeira, “Anisotropic voter model,” Journal of statistical physics, vol. 78, no. 3-4, pp. 963–970, 1995.
  * (19) A. R. Vieira and N. Crokidakis, “Phase transitions in the majority-vote model with two types of noises,” Physica A: Statistical Mechanics and its Applications, vol. 450, pp. 30–36, 2016.
  * (20) A. L. M. Vilela and F. G. B. Moreira, “Majority-vote model with different agents,” Physica A: Statistical Mechanics and its Applications, vol. 388, p. 4171, 2009.
  * (21) A. L. M. Vilela and H. E. Stanley, “Effect of strong opinions on the dynamics of the majority-vote model,” Scientific reports, vol. 8, no. 1, pp. 1–8, 2018.
  * (22) A. L. M. Vilela, L. F. C. Pereira, L. Dias, H. E. Stanley, and L. R. da Silva, “Majority-vote model with limited visibility: An investigation into filter bubbles,” Physica A: Statistical Mechanics and its Applications, vol. 563, p. 125450, 2021.
  * (23) J. H. Feldhoff, S. Lange, J. Volkholz, J. F. Donges, J. K. Jürgen, and F.-W. Gerstengarbe, “Complex networks for climate model evaluation with application to statistical versus dynamical modeling of south american climate,” Climate Dynamics, vol. 44, no. 5-6, pp. 1567–1581, 2015.
  * (24) J. C. Reijneveld, S. C. Ponten, H. W. Berendse, and C. J. Stam, “The application of graph theoretical analysis to complex networks in the brain,” Clinical neurophysiology, vol. 118, no. 11, pp. 2317–2331, 2007.
  * (25) A.-L. Barabási, R. Albert, and H. Jeong, “Scale-free characteristics of random networks: the topology of the world-wide web,” Physica A: Statistical Mechanics and its Applications, vol. 281, no. 1-4, pp. 69–77, 2000\.
  * (26) X. lei An, L. Zhang, Y. zhen Li, and J. gang Zhang, “Synchronization analysis of complex networks with multi-weights and its application in public traffic network,” Physica A: Statistical Mechanics and its Applications, vol. 412, pp. 149–156, 2014.
  * (27) P. Erdös, A. Rényi, et al., “On the evolution of random graphs,” Publ. Math. Inst. Hung. Acad. Sci, vol. 5, no. 1, pp. 17–60, 1960\.
  * (28) B. Bollobás and A. Thomason, “Random graphs of small order,” in North-Holland Mathematics Studies, vol. 118, pp. 47–97, Elsevier, 1985.
  * (29) R. Cont, J.-P. Bouchaud, et al., “Herd behavior and aggregate fluctuations in financial markets,” Macroeconomic Dynamics, vol. 4, no. 2, pp. 170–196, 2000.
  * (30) K. Sznajd-Weron and R. Weron, “A simple model of price formation,” International Journal of Modern Physics C, vol. 13, no. 01, pp. 115–123, 2002\.
  * (31) T. Lux and M. Marchesi, “Scaling and criticality in a stochastic multi-agent model of a financial market,” Nature, vol. 397, no. 6719, pp. 498–500, 1999.
  * (32) S. Bornholdt, “Expectation bubbles in a spin model of markets: Intermittency from frustration across scales,” International Journal of Modern Physics C, vol. 12, no. 05, pp. 667–674, 2001.
  * (33) T. Kaizoji, S. Bornholdt, and Y. Fujiwara, “Dynamics of price and trading volume in a spin model of stock markets with heterogeneous agents,” Physica A: Statistical Mechanics and its Applications, vol. 316, no. 1-4, pp. 441–452, 2002.
  * (34) T. Takaishi, “Simulations of financial markets in a potts-like model,” International Journal of Modern Physics C, vol. 16, no. 08, pp. 1311–1317, 2005\.
  * (35) B. Mandelbrot, “The Variation of Certain Speculative Prices,” The Journal of Business, vol. 36, pp. 394–394, 1963.
  * (36) M. Mazura, M. Dang, and M. Vega, “Covid-19 and the march 2020 stock market crash. evidence from s&p1500,” Finance Research Letters, vol. 38, p. 101690, 2021.
  * (37) S. Fu, C. Liu, and X. Wei, “Contagion in global stock markets during the covid-19 crisis,” Global Challenges, vol. 5, no. 10, p. 2000130, 2021.


# The effect of disorder on phases across two-dimensional thermal melting

Prashanti Jami  Department of Physics, Indian Institute of Science Education
and Research (IISER) Kolkata, Mohanpur - 741246, West Bengal, India  Pinaki
Chaudhuri  The Institute of Mathematical Sciences, Taramani, Chennai 600113,
India  Chandan Dasgupta  Department of Physics, Indian Institute of Science,
Bangalore 560012, India  International Centre for Theoretical Sciences, Tata
Institute of Fundamental Research, Bangalore 560089, India  Amit Ghosal
Department of Physics, Indian Institute of Science Education and Research
(IISER) Kolkata, Mohanpur - 741246, West Bengal, India

###### Abstract

We study melting in a two-dimensional system of classical particles with
Gaussian-core interactions in disordered environments. The pure system
validates the conventional two-step melting with a hexatic phase intervening
between the solid and the liquid. This picture is modified in the presence of
pinning impurities. A random distribution of pinning centers forces a hexatic-
like low temperature phase that transits into a liquid at a single melting
temperature $T^{\rm RP}_{\rm m}$. In contrast, pinning centers located at
randomly chosen sites of a perfect crystal anchors a solid at low temperatures
which undergoes a direct transition to the liquid at $T^{\rm CP}_{\rm m}$.
Thus, the two-step melting is lost in either cases of disorder. We discuss the
characteristics of melting depending on the nature of the impurities.

Introduction — Enhanced fluctuations make two-dimensional melting a topic of
immense research interest. Unlike their three-dimensional counterparts
undergoing “Lindemann melting” [1, 2], 2D melting is mediated by the unbinding
of topological defects. The positional order (PO) and bond-orientational order
(BOO) decouple in 2D, producing a “hexatic phase” sandwiched between the solid
and the liquid. Hexaticity, a rich concept, is realized in colloids [3, 4],
the vortex lattice in superconductors [5], in active Brownian disks [6] and
recently in van der Waals magnet [7]. The celebrated KTHNY theory [8, 9, 10,
11, 12], pictures 2D-melting as a two-step process involving successive
unbinding of dislocations and disclinations– presented schematically in Fig.
1(a). However, the relevance of the two-step 2D melting has also been debated
[13, 14, 15, 16].

Figure 1: A schematic representation of 2D-melting in: (a) KTHNY description
in a pure system, (b) A system with randomly placed impurities, and (c) With a
fraction of particles frozen at randomly chosen zero-temperature positions of
the pure system, suggesting that the clean-melting may be obscured by
impurities (See text).

Quenched disorder, inherent to real materials, can not only move around the
phase boundaries but is also capable of modifying the mechanism of melting.
For example, impurities can generate unbounded defects even at $T=0$, and
thereby mask the unbinding of thermal defect-pairs. This could strike out
solidity even at the lowest $T$, as suggested by Nelson [17] – portrayed
schematically in Fig. 1(b). In contrast, impurities which pin a given fraction
of particles on sites of the underlying perfect lattice could stabilize the
solid by anchoring it via these commensurate locations, and thereby consume
the phase space of hexaticity (see Fig. 1(c)). Role of disorder in
destabilizing the hexatic phase [18] and in enhancing long-range correlation
has also been pointed out [19]. Thus, a careful analysis of 2D-melting in
disordered media can potentially uncover new paradigms.

Experiments on colloids [20], vortex lattices [21, 22] and multi component
mixtures [23] indicate a broadened stability of the hexatic phase in the
presence of disorder consistent with recent calculations [24, 25, 26, 27, 28,
32]. Study of 2D melting is also popular in confined geometry [29, 30, 31]
mimicking disordered background. Zeng et al.[33] have argued that a solid
(“Bragg glass”) phase with power-law decay of translational correlations
cannot occur in a 2D system with impurities. Pronounced hexatic correlations
are expected to be present [34, 35] if the disorder is not strong, though
there are controversies [36, 37, 38] about the existence of a hexatic glass
phase with long- or quasi-long-range hexatic order in 2D. In contrast, phase
boundaries of the two-step melting are found to be insensitive to quenched
defects on a spherical surface [39].

In this letter, we investigate the phases across melting of a bulk 2D system
of soft-core particles, modeled via Gaussian-core interactions [40], which is
known to validate the KTHNY-melting scenario in a pure system [41]. Addressing
the role of quenched disorder in the phase behaviour this model, our key
results are summarized as follows: (i) Random-pinning (RP) destabilizes
solidity causing a single transition from a low-$T$ hexatic-like phase to a
high-$T$ liquid. Here, the low-$T$ phase undergoes a likely crossover from
hexatic-glass to hexatic-liquid. (ii) On the other hand, the commensurate-
pinning (CP) anchors solidity and engulfs hexaticity – even the high-$T$
liquid phase supports inhomogeneous pockets of crystallinity. The defect
locations correlate oppositely with pinning centers in the two models of
disorder – defects tend to bind with the pinning centers for RP-systems,
whereas they stay away from the impurities in CP-systems. Thus, in either
realization of the quenched disorder, the two-step melting is lost.

Model & method: — We introduce disorder in two different ways: (a) Random
pinning (RP), in which we freeze a given fraction ($n_{\rm imp}$) of
particles, chosen randomly in space, within a high-$T$ liquid configuration.
Here, these immobile particles act as disorder. (b) Commensurate pinning (CP),
where $n_{\rm imp}$ fraction of particles are frozen at randomly chosen
positions of an ideal triangular lattice – the ground state configuration of
the pure system. Note that CP represents correlated disorder with a long-range
positional correlation of a perfect lattice. In contrast, RP constitutes
nearly uncorrelated disorder though weak short-range correlation of a high-$T$
liquid may exist. We investigated a system of $N=4356$ particles, with $n_{\rm
imp}=3.5\%$. These results were compared with those from a pure system with
$N=4096$ particles. For these systems, we sample configurations via molecular
dynamics. [42] within the canonical ensemble, using LAMMPS [43] We consider a
simulation box having dimensions $L_{x}=\frac{2}{\sqrt{3}}L_{y}$, having
periodic boundary conditions. $L_{x}$ is adjusted to keep the density, $\rho$
of particles fixed for all our studies ($\rho=0.628$). We carried out $2\times
10^{7}$ MD steps with a time step $\delta t=0.005$. We use dimensionless
parameters: $t^{\prime}=t\sqrt{\epsilon/m\sigma^{2}}$ and
$E^{\prime}=E/\epsilon$, where $m$ is the mass of each particle. $T$ is
expressed in terms of $\Gamma^{-1}$[44], where
$\Gamma=\epsilon\exp(-\sqrt{3}/2\rho)/K_{B}T$. The physical observables are
averaged over $8$ \- $10$ independent pinning configurations for a given
$n_{\rm imp}$.

Positional and bond orientational order: — A pure 2D solid is characterized by
two kinds of ordering: (i) PO measured by $\psi_{\rm
T}=\frac{1}{N}\left\langle|\Psi_{\rm T}|\right\rangle$, where $\Psi_{\rm
T}=\sum_{i=1}^{N}\exp(i{\bf G}.{\bf r}_{i})$. $\bf{G}$ is a first shell
reciprocal-lattice vector of the underlying triangular crystal and ${\bf
r}_{i}$ is the position of particle $i$, and (ii) BOO, quantified by
$\psi_{\rm 6}=\frac{1}{N}\left\langle|\sum_{k=1}^{N}\Psi_{\rm
6}(r_{k})|\right\rangle$ where $\Psi_{\rm
6}=\frac{1}{N_{b}(k)}\sum_{l=1}^{N_{b}(k)}\exp{(i6\theta_{kl})}$. The sum is
over the $N_{b}(k)$ nearest neighbors of particle $k$ identified by a Voronoi
construction [45] and $\theta_{kl}$ is the angle that a line joining particle
$k$ and particle $l$ makes with a reference axis. KTHNY theory predicts two
critical temperatures: $\Gamma^{-1}_{\rm SH}$ and $\Gamma^{-1}_{\rm HL}$, for
the thermal depletion of quasi-long-range PO (solid to hexatic) and BOO
(hexatic to liquid) respectively, leaving a hexatic phase with quasi-long-
range BOO between the solid and isotropic liquid phases.

In Fig. 2(a,b), we plot the thermal evolution of $\psi_{\rm T}$ and $\psi_{\rm
6}$ for pure, RP and CP systems. While the pure system follows KTHNY melting
111The three phases are clearly identified by the snapshots at three
representatives $T$s in the SM, Fig. S1(a-c) with $\Gamma^{-1}_{\rm
SH}=0.0140$ and $\Gamma^{-1}_{\rm HL}=0.0162$, $\psi_{\rm T}$ in the CP-system
is found to survive to larger $T$. The RP-system shows a much weaker
$\psi_{\rm T}$ than the other two, even at the lowest $T$ and depletes very
gradually with $T$ without any threshold behavior. A threshold behavior near
$\Gamma^{-1}_{\rm HL}$ is also seen in the pinned systems in Fig. 2(b), albeit
the transitions are broader. Unlike $\psi_{\rm T}$, the $\psi_{\rm 6}$ is
comparable at low-$T$ in pure, CP- and RP-systems. We also note that
$\psi_{\rm T}$ and $\psi_{\rm 6}$ show a significant drop at the same critical
temperature in a CP-system, implying a direct transition from solid to liquid,
which we discuss further below.

In addition, the fluctuations of $\psi_{\rm T}$ and $\psi_{\rm 6}$ define
generalised susceptibilities
$\chi_{\alpha}=\frac{1}{N}\left[\langle|\Psi_{\alpha}|^{2}\rangle-\langle|\Psi_{\alpha}|\rangle^{2}\right]$,
(with $\alpha=T$ or $6$), and help to identify $\Gamma^{-1}_{\rm SH}$ and
$\Gamma^{-1}_{\rm HL}$, as shown in Fig. 2(c,d). Their behavior confirms that
the pure system shows sharp transitions. Consistent with our finding in panel
(a), $\chi_{\rm T}$ in the RP-system features only a broad and low hump,
hinting that a low-$T$ phase in such a system represents a broad crossover
between a hexatic glass [34] and a hexatic liquid 222The low-$T$ phase support
significant hexatic order $\psi_{\rm 6}$ for our model parameters, while
$\psi_{\rm T}$ and the snapshots establish the amorphous and glassy nature.
The $\Gamma^{-1}$-dependence of $\psi_{\rm 6}$ and $\psi_{\rm T}$ is
indicative of a crossover from a hexatic-glass to a hexatic-liquid, before the
RP-system transits to a liquid at $\Gamma^{-1}_{\rm RP}$. However, the
resolution of our simulation is inadequate for drawing a firm conclusion..
This is also consistent with the trajectory picture of RP system at
$\Gamma^{-1}=0.0028$ in the supplementary material (SM) [48] Fig. S1(d).
Congruous with our findings in panels (a,b), the locations of the peak of
$\chi_{\rm T}$ and $\chi_{\rm 6}$ verify that PO and BOO in the CP-system
vanish at a single $\Gamma^{-1}_{\rm CP}$. While our results from Fig. 2 seem
to support the schematic phase-diagram of Fig. 1, we emphasize that the
‘impure’ phases at low- and high-$T$ defy conventional wisdom. These include
the presence of unbound defects even at $T=0$ in RP-systems, and pockets of
crystallinity deep into the liquid phase in CP-systems, as seen from Fig.
S1(g-i) in the SM [48], and discussed below.

Figure 2:  Positional and orientational ordering tendencies are shown as a
function of $\Gamma^{-1}$ for pure (green circles), RP (red triangles) and CP
(blue square) systems. (a) Decay of PO ($\psi_{\rm T}$) with $\Gamma^{-1}$
(the inset shows $\psi_{\rm T}$ for $\Gamma^{-1}\rightarrow 0$). (b) The
softening of the BOO ($\psi_{\rm 6}$). Here $\psi_{\rm
6}(\Gamma^{-1}\rightarrow 0)$ is shown in inset. Panel (c) and (d) show the
corresponding susceptibilities $\chi_{\rm T}$ and $\chi_{\rm 6}$. The location
of the peaks in $\chi_{\rm T}$ and $\chi_{\rm 6}$ identify the transitions. We
find $\Gamma^{-1}_{\rm RP}=0.0145$ from $\chi_{\rm 6}$. Interestingly, the
peak in $\chi_{\rm T}$ and $\chi_{\rm 6}$ for CP-systems appear at the same
$\Gamma^{-1}_{\rm CP}=0.0187$.  Figure 3:  The distribution of distances
between dislocation pairs, $P(r_{\rm dp})$ in different phases. Panels (a)-(c)
show results for pure, RP- and CP-systems, respectively, at low-$T$. $a_{0}$
is the average inter-particle distance. $P(r_{\rm dp})$ is sharply peaked at
the typical distance between dislocation pairs, though there are differences
in details. The distribution has a short tail for the pure system, a very long
tail for RP-systems, and essentially no tail for CP-systems. Panel (d)
displays $\Delta P(r_{\rm dp})$ taken in pure system between $T^{\prime}s$
just above and below $\Gamma^{-1}_{\rm SH}$. The inset shows corresponding
$P(r_{\rm dp})$. The largest zero-crossing distance of $\Delta P(r_{\rm dp})$
is taken as $r^{c}_{\rm dp}$, and is marked as an arrow in panel (a). Panels
(e) and (f) show the spatial density of unbound defects just prior to melting
into a liquid for a given realization of RP- and CP-system. The pinning
centers are marked via magenta dots. Figure 4: Evolution of defects: Panels
(a)-(c) represent the variation of defects with $\Gamma^{-1}$ (green: pure,
red: RP, and blue: CP-system). (a) The evolution of unbound dislocations (in
squares) as well as disclinations (in circles) validate KTHNY melting in pure
systems; (traces of disclinations are multiplied by $10$ for visual clarity).
(b) The presence of significant unbound dislocations in RP-systems at low
$\Gamma^{-1}$ prohibits solidity. In CP-systems, shown in panel (c), the
proliferation of unbound defects (both dislocations and disclinations)
commences at a single threshold ($\Gamma^{-1}$=0.0198), implying a direct
solid-liquid transition. Panel (d) shows the comparative number defects in RP-
systems with reference to a pure system, illustrating how the number of bound
(B) and unbound (UB) defects grow with increasing temperature.

Defects analysis: — In KTHNY theory, a 2D-solid transits to the hexatic phase
by the unbinding of paired dislocations [10, 11]. We proceed to examine the
consistency of this picture in Fig. 2. This requires an estimation of the
critical distance between two dislocations (with equal and opposite Burger
vectors) below which they are bound. We first employ the Hungarian algorithm
[49], which chooses ‘correct’ partners of dislocations by minimizing the sum
of the distances between all partners [50]. The distribution of resulting pair
distances, $P(r_{\rm dp})$, at low-$T$ is presented in Fig. 3(a-c). $P(r_{\rm
dp})$ is sharply peaked for pure systems in Fig. 3(a), with insignificant
weight at larger $r_{\rm dp}$. In contrast, its long tail at low $T$ (shown
for $\Gamma^{-1}=0.0056$ in Fig. 3(b)) for the RP-system arises from unbound
dislocations even for $T\rightarrow 0$, which destabilize a true solid.
$P(r_{\rm dp})$ in CP-systems (Fig. 3(c)) consists of the initial sharp peak,
and nearly no weight for larger $r_{\rm dp}$. A discernible tail in $P(r_{\rm
dp})$ for pure systems develops when dislocation pairs start unbinding. An
integrated distribution of $P(r_{\rm dp})$ features a threshold behavior at
this transition (see SM [48]).
To obtain the critical $r_{\rm dp}$ for the pure system, we plot in Fig. 3(d)
the difference of these distributions, $\Delta P(r_{\rm dp})$, at temperatures
just above and below $\Gamma^{-1}_{\rm SH}$, while the corresponding $P(r_{\rm
dp})$-s are shown as the inset. The total positive and negative weights of
$\Delta P(r_{\rm dp})$ cancel out, and $r^{c}_{\rm dp}\approx 2.15a_{0}$ is
identified as the last zero-crossing point. This identification is found
robust for $T$’s near $\Gamma^{-1}_{\rm SH}$. A study of distances of
disclination pairs yielded a similar critical distance between disclination
pairs. Once extracted for the pure system, these critical distances were used
for analyzing pinned systems. Subsequently, we explored the thermal evolution
of the defects and their unbinding in Fig. 4(a-c). For the lowest $T$, defects
are essentially absent in the pure system. Unbound disclinations proliferate
at $\Gamma^{-1}\approx 0.0162$, whereas dislocations unbind at
$\Gamma^{-1}\approx 0.0145$, with a hexatic phase at intervening temperatures
[41], consistent with Fig. 2. The CP-system (Fig. 4(c)) behaves like a
‘better’ solid at low-$T$ due to the absence of any free defects up to
$\Gamma^{-1}=0.0182$ beyond that unbound dislocations and disclinations start
proliferating at the same $\Gamma^{-1}_{\rm CP}$. There are significant number
of impurity-induced unbound dislocations in the RP-system for $T\rightarrow
0$, as also concluded from Fig. 2. Here, unpaired dislocations are not only
present for all $T$, they even outnumber bound dislocations at low-$T$. Fig.
4(d) addresses the role of the impurity induced free defects (at $T\rightarrow
0$) in RP-systems, on the thermal defects, whose unbinding drives the two
transitions in a pure system. Number of bound (B) and unbound (UB) defects,
with corresponding numbers subtracted for an equivalent pure system, are
examined separately in Fig. 4(d). These numbers increase sharply with $T$
until the system transits to the liquid. Thus, the impurity induced defects
help promote further thermal defects, than in pure systems for
$\Gamma^{-1}<\Gamma^{-1}_{\rm HL}$. Such a rise disappears in the liquid. In
fact, this difference in bound defects in the liquid goes down to a even lower
value than the corresponding number at $T\rightarrow 0$.

Figure 5: Correlations: Panels (a), (b) show the orientational correlation
function, $g_{6}(r)$, of RP- and CP-systems, respectively, for the values of
$\Gamma^{-1}$ as labelled in each case. In panel (a), power law decay of the
correlations is visible (for $\Gamma^{-1}<{}^{\rm RP}\Gamma^{-1}_{\rm m}$)
down to a very low $T$ in the RP-system; e.g. we find
$g_{6}(r,\Gamma^{-1}=0.0112)\sim(r/a_{0})^{-0.14}$,
$g_{6}(r,\Gamma^{-1}=0.0.014)\sim(r/a_{0})^{-0.2}$,
$g_{6}(r,\Gamma^{-1}=0.0154)\sim(r/a_{0})^{-0.64}$, while
$g_{6}(r,\Gamma^{-1}=0.0162)\sim{\rm exp}^{-0.16(r/a_{0})}$, as shown via
dashed lines. In contrast, in CP-system, as shown in (b), correlations get
modified. Even in the liquid phase, the conventional exponential decay
flattens out at large $r$, implying a “remnant crystallinity” (See text).

Correlations: — Finally, we discuss the orientational correlations as measured
by the correlation function $g_{6}(r)=\langle\Psi_{6}({\bf r}_{i})\Psi_{\rm
6}^{*}({\bf r}_{j})\rangle$ where $r=\lvert{\bf r}_{i}-{\bf r}_{j}\rvert$. The
$T$-dependence of orientational correlations in the pure system follows KTHNY
scenario [10, 11, 12] as claimed earlier [41]. The evolution of $g_{6}(r)$ for
various $T$ is shown in Fig. 5 for RP- and CP-systems. Our
$\chi^{2}$-minimization analysis [51] of the large-$r$ decay of $g_{6}(r)$ in
the RP-system (Fig. 5(a)) identified a power-law behavior for nearly the
entire low-$T$ phase. This power-law behavior continues until an exponential
decay sets in for $\Gamma^{-1}\geq 0.0162$, signaling the onset of liquidity.
Intriguingly, $g_{6}(r)$ in CP-systems shows the enhanced solidity for
$\Gamma^{-1}\leq 0.0187$, where its traces remain largely flat. Beyond the
direct melting from solid to liquid for $\Gamma^{-1}>0.0190$, $g_{6}(r)$ in
CP-systems displays a tendency of plateauing at large-$r$, though it decays at
intermediate $r$. This is a signature of ‘remnant solidity’ arising from local
crystalline pockets surrounding the impurities whose locations are
commensurate with the perfect crystal and hence anchoring crystallinity in the
vicinity (See SM [48]). This is a direct consequence of the correlated nature
of the CP-impurities.
Conclusion: — To summarize, we demonstrate that the conventional picture of
2D-melting undergoes significant changes in the presence of impurities. While
RP-disorder destabilizes solidity and CP-disorder removes the hexatic phase,
the low-$T$ phase in RP-systems is not the conventional hexatic. Similarly,
the high-$T$ phase in the CP-systems mixes remnant solidity with the liquid
phase. The inhomogeneous melting (Fig. S1 in SM [48] generates defects which
correlate differently with pinning centers: For RP-systems the defects tend to
bind with the pinning centers, whereas the defects stay away from the
impurities in CP-systems. While defects are found essential for driving the
melting, our MD configurations indicate that they often bunch up in various
shapes of macroscopic size (See videos in SM [48]). An extension of our study
to larger systems exploring possible role of grain boundaries on melting is a
promising future direction. It will also be interesting to inspect the role of
quantum fluctuations in these thermal phases. We hope that our findings will
motivate future experiments for shedding new light.

## References

  * Lindemann [1910] F. Lindemann, The calculation of molecular vibration frequencies phys, Z 11, 609 (1910).
  * Lozovik [1987] Y. E. Lozovik, Ion and electron clusters, Soviet Physics Uspekhi 30, 912 (1987).
  * Zahn  _et al._ [1999] K. Zahn, R. Lenke, and G. Maret, Two-stage melting of paramagnetic colloidal crystals in two dimensions, Phys. Rev. Lett. 82, 2721 (1999).
  * Gasser  _et al._ [2010] U. Gasser, C. Eisenmann, G. Maret, and P. Keim, Melting of crystals in two dimensions, ChemPhysChem 11, 963 (2010).
  * Guillamón  _et al._ [2009] I. Guillamón, H. Suderow, A. Fernández-Pacheco, J. Sesé, R. Córdoba, J. M. De Teresa, M. R. Ibarra, and S. Vieira, Direct observation of melting in a two-dimensional superconducting vortex lattice, Nature Physics 5, 651 (2009).
  * Digregorio  _et al._ [2018] P. Digregorio, D. Levis, A. Suma, L. F. Cugliandolo, G. Gonnella, and I. Pagonabarraga, Full phase diagram of active brownian disks: From melting to motility-induced phase separation, Phys. Rev. Lett. 121, 098003 (2018).
  * Meisenheimer  _et al._ [2023] P. Meisenheimer, H. Zhang, D. Raftrey, X. Chen, Y.-T. Shao, Y.-T. Chan, R. Yalisove, R. Chen, J. Yao, M. C. Scott, W. Wu, D. A. Muller, P. Fischer, R. J. Birgeneau, and R. Ramesh, Ordering of room-temperature magnetic skyrmions in a polar van der waals magnet, Nature Communications 14, 3744 (2023).
  * Kosterlitz and Thouless [1972] J. M. Kosterlitz and D. J. Thouless, Long range order and metastability in two dimensional solids and superfluids. (application of dislocation theory), Journal of Physics C: Solid State Physics 5, L124 (1972).
  * Kosterlitz and Thouless [1973] J. M. Kosterlitz and D. J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems, J. Phys. C 6, 1181 (1973).
  * Halperin and Nelson [1978] B. I. Halperin and D. R. Nelson, Theory of two-dimensional melting, Phys. Rev. Lett. 41, 121 (1978).
  * Nelson and Halperin [1979] D. R. Nelson and B. I. Halperin, Dislocation-mediated melting in two dimensions, Phys. Rev. B 19, 2457 (1979).
  * Young [1979] A. P. Young, Melting and the vector coulomb gas in two dimensions, Phys. Rev. B 19, 1855 (1979).
  * Kapfer and Krauth [2015] S. C. Kapfer and W. Krauth, Two-dimensional melting: From liquid-hexatic coexistence to continuous transitions, Phys. Rev. Lett. 114, 035702 (2015).
  * Chui [1983] S. T. Chui, Grain-boundary theory of melting in two dimensions, Phys. Rev. B 28, 178 (1983).
  * Qi  _et al._ [2014] W. Qi, A. P. Gantapara, and M. Dijkstra, Two-stage melting induced by dislocations and grain boundaries in monolayers of hard spheres, Soft Matter 10, 5449 (2014).
  * Mazars [2015] M. Mazars, The melting of the classical two-dimensional wigner crystal, Europhysics Letters 110, 26003 (2015).
  * Nelson [1983] D. R. Nelson, Reentrant melting in solid films with quenched random impurities, Phys. Rev. B 27, 2902 (1983).
  * Qi and Dijkstra [2015] W. Qi and M. Dijkstra, Destabilisation of the hexatic phase in systems of hard disks by quenched disorder due to pinning on a lattice, Soft Matter 11, 2852 (2015).
  * Guillamón  _et al._ [2014] I. Guillamón, R. Córdoba, J. Sesé, J. M. De Teresa, M. R. Ibarra, S. Vieira, and H. Suderow, Enhancement of long-range correlations in a 2d vortex lattice by an incommensurate 1d disorder potential, Nature Physics 10, 851 (2014).
  * Deutschländer  _et al._ [2013] S. Deutschländer, T. Horn, H. Löwen, G. Maret, and P. Keim, Two-dimensional melting under quenched disorder, Phys. Rev. Lett. 111, 098301 (2013).
  * Ganguli  _et al._ [2016] S. C. Ganguli, H. Singh, I. Roy, V. Bagwe, D. Bala, A. Thamizhavel, and P. Raychaudhuri, Disorder-induced two-step melting of vortex matter in co-intercalated $\mathrm{NbS}{\mathrm{e}}_{2}$ single crystals, Phys. Rev. B 93, 144503 (2016).
  * Duhan  _et al._ [2023] R. Duhan, S. Sengupta, R. Tomar, S. Basistha, V. Bagwe, C. Dasgupta, and P. Raychaudhuri, Structure and dynamics of a pinned vortex liquid in superconducting $a$-rexzr (x $\approx$ 6) thin film (2023), arXiv:2304.10926 [cond-mat.supr-con] .
  * Li  _et al._ [2023] Y.-W. Li, Y. Yao, and M. P. Ciamarra, Two-dimensional melting of two- and three-component mixtures, Phys. Rev. Lett. 130, 258202 (2023).
  * Tsiok  _et al._ [2021] E. N. Tsiok, Y. D. Fomin, E. A. Gaiduk, and V. N. Ryzhov, Structural transition in two-dimensional hertzian spheres in the presence of random pinning, Phys. Rev. E 103, 062612 (2021).
  * Gaiduk  _et al._ [2019] E. A. Gaiduk, Y. Fomin, E. N. Tsiok, and V. N. Ryzhov, The influence of random pinning on the melting scenario of two-dimensional soft-disk systems, Molecular Physics 117, 2910 (2019), https://doi.org/10.1080/00268976.2019.1607917 .
  * Shankaraiah  _et al._ [2020] N. Shankaraiah, S. Sengupta, and G. I. Menon, Disorder-induced enhancement of local hexatic correlations in two-dimensional fluids, Journal of Physics: Condensed Matter 32, 184003 (2020).
  * Arjun H and Chaudhuri [2020] Arjun H and P. Chaudhuri, Dense hard disk ordering: influence of bidispersity and quenched disorder, Journal of Physics: Condensed Matter 32, 414001 (2020).
  * Tsiok  _et al._ [2015] E. N. Tsiok, D. E. Dudalov, Y. D. Fomin, and V. N. Ryzhov, Random pinning changes the melting scenario of a two-dimensional core-softened potential system, Phys. Rev. E 92, 032110 (2015).
  * Melzer  _et al._ [2012] A. Melzer, A. Schella, T. Miksch, J. Schablinkski, D. Block, A. Piel, H. Thomsen, H. Kählert, and M. Bonitz, Phase transitions of finite dust clusters in dusty plasmas, Contributions to Plasma Physics 52, 795 (2012), https://onlinelibrary.wiley.com/doi/pdf/10.1002/ctpp.201200028 .
  * Ash  _et al._ [2017] B. Ash, J. Chakrabarti, and A. Ghosal, Static and dynamic properties of two-dimensional coulomb clusters, Phys. Rev. E 96, 042105 (2017).
  * Ash  _et al._ [2018] B. Ash, C. Dasgupta, and A. Ghosal, Analysis of vibrational normal modes for coulomb clusters, Phys. Rev. E 98, 042134 (2018).
  * Cha and Fertig [1995] M.-C. Cha and H. A. Fertig, Disorder-induced phase transitions in two-dimensional crystals, Phys. Rev. Lett. 74, 4867 (1995).
  * Zeng  _et al._ [1999] C. Zeng, P. L. Leath, and D. S. Fisher, Absence of two-dimensional bragg glasses, Phys. Rev. Lett. 82, 1935 (1999).
  * Chudnovsky [1989] E. M. Chudnovsky, Hexatic vortex glass in disordered superconductors, Phys. Rev. B 40, 11355 (1989).
  * Chudnovsky [1991a] E. M. Chudnovsky, Orientational and positional order in flux lattices of type-ii superconductors, Phys. Rev. B 43, 7831 (1991a).
  * Toner [1991a] J. Toner, Orientational order in disordered superconductors, Phys. Rev. Lett. 66, 2523 (1991a).
  * Chudnovsky [1991b] E. M. Chudnovsky, Comment on “orientational order in disordered superconductors.”, Phys. Rev. Lett. 67, 1809 (1991b).
  * Toner [1991b] J. Toner, Toner replies, Phys. Rev. Lett. 67, 1810 (1991b).
  * Singh  _et al._ [2022] N. Singh, A. Sood, and R. Ganapathy, Observation of two-step melting on a sphere, Proceedings of the National Academy of Sciences 119, e2206470119 (2022).
  * Stillinger [2008] F. H. Stillinger, Phase transitions in the Gaussian core system, The Journal of Chemical Physics 65, 3968 (2008), https://pubs.aip.org/aip/jcp/article-pdf/65/10/3968/11360037/3968_1_online.pdf .
  * Prestipino  _et al._ [2011] S. Prestipino, F. Saija, and P. V. Giaquinta, Hexatic phase in the two-dimensional gaussian-core model, Phys. Rev. Lett. 106, 235701 (2011).
  * Frenkel and Smit [2001] D. Frenkel and B. Smit,  _Understanding molecular simulation: from algorithms to applications_ , Vol. 1 (Elsevier, 2001).
  * Plimpton [1995] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, Journal of Computational Physics 117, 1 (1995).
  * Gann  _et al._ [1979] R. C. Gann, S. Chakravarty, and G. V. Chester, Monte carlo simulation of the classical two-dimensional one-component plasma, Phys. Rev. B 20, 326 (1979).
  * Tipper [1991] J. C. Tipper, Fortran programs to construct the planar voronoi diagram, Computers & Geosciences 17, 597 (1991).
  * Note [1] The three phases are clearly identified by the snapshots at three representatives $T$s in the SM, Fig. S1(a-c).
  * Note [2] The low-$T$ phase support significant hexatic order $\psi_{\rm 6}$ for our model parameters, while $\psi_{\rm T}$ and the snapshots establish the amorphous and glassy nature. The $\Gamma^{-1}$-dependence of $\psi_{\rm 6}$ and $\psi_{\rm T}$ is indicative of a crossover from a hexatic-glass to a hexatic-liquid, before the RP-system transits to a liquid at ${}^{\rm RP}\Gamma^{-1}_{\rm m}$. However, the resolution of our simulation is inadequate for drawing a firm conclusion.
  * [48] The URL for SM will appear here.
  * Kuhn [1955] H. W. Kuhn, The hungarian method for the assignment problem, Naval research logistics quarterly 2, 83 (1955).
  * Lau and Dasgupta [1989] M.-h. Lau and C. Dasgupta, Numerical investigation of the role of topological defects in the three-dimensional heisenberg transition, Phys. Rev. B 39, 7212 (1989).
  * Levenberg [1944] K. Levenberg, A method for the solution of certain non-linear problems in least squares, Quarterly of applied mathematics 2, 164 (1944).

Supplementary material for ‘The effect of disorder on phases across two-
dimensional thermal melting’

In order to support the key conclusions reported in the main manuscript, we
include below additional results in this supplementary materials (SM).

##  I Model and Methods

For our study, we consider the Gaussian-core model in a 2D system with
periodic boundary conditions, within which particles act via the following
interaction. $V(r)=\epsilon exp(-r^{2}/\sigma^{2})$, with $\epsilon>0$,
$\epsilon$ and $\sigma$ are the energy and length scales, respectively. We
carry out molecular dynamics (MD) simulations of this model system, using
LAMMPS, starting with a high temperature and cooling the system down to the
desired temperature, as detailed in the main text. We use $10^{5}$ MD steps
($t=5000$) with sampling time window of $t=0.1$, which are recorded after
$10^{7}$ sweeps of equilibration runs. The desired temperatures are maintained
via the Berendsen thermostat. To ensure correct equilibration before taking
statistics, we looked into the distribution of velocities of the particles,
which assumes the form of an Maxwell-Boltzmann distribution in thermal
equilibrium as well as independence of the temporal correlation of the
observables on the time origin, i.e.
$\zeta(t_{1},t_{2})=\zeta_{(}t_{1}-t_{2})$, where $\zeta$ is any temporal
correlation defined between two time points $t_{1}$ and $t_{2}$.

Figure S1: Trajectories: The trajectories of the particles studied at
different temperatures ($\Gamma^{-1}$) for all three systems (the pinning
centers marked in red) for $\Delta t=5000$. Panels (a-c) represent the
trajctories of particles in the clean system. For lowest $\Gamma^{-1}$ (panel
a), the particles remain firmly tied to the equilibrium position of a
triangular lattice reflecting the solid phase. In the hexatic phase at
$\Gamma^{-1}=0.0145$ (panel b), as per Fig. 2 of the main text, particles tend
to move preferentially along the principal direction. In the liquid phase
(panel c), particles execute isotropic motion, covering the entire space. The
RP-system at low-$T$ (panel d) features distorted lattice lines and a
resulting irregular array of particles. Together with the results from Fig. 2
of main text, the trajectory of particles indicate that the lowest-$T$ phase
of RP-system constitute a hexatic glass. Panel (e) depicts the trajectory plot
near ${}^{\rm RP}\Gamma^{-1}_{\rm m}$ and features incipient melting which
results into motional signatures which are rather inhomogeneous in space.
Finally, the liquid phase in the RP-system in panel (f) displays particles’
motion everywhere except for a small region surrounding the impurities. The
motion in CP-system looks similar to the clean system in the solid phase
(panel g). Close to the ${}^{\rm CP}\Gamma^{-1}_{\rm m}$ in panel (h) the
dynamics shows inhomogeneity, but unlike in RP-system the mobile particles
tend to stay away from pinning centers. the concentration of pinned centers
keeps away the melting from the pinned centers (panel h). While most particles
delocalize significantly in the high-$T$ phase of the CP-system (panel i),
local patches of solid-like regions still survive, as seen from the trajectory
plot.

##  II Snapshots of particle trajectories

Here we plot the trajectories of particles during their equilibrium dynamics
for the clean, RP- and CP-systems. Such motional signatures are displayed in
Fig.S1 for three representative temperatures, one at low-$T$ in the solid
phase (likely a hexatic glass in case of RP-system), at intermediate
temperature (just before the system melts into an isotropic liquid (in hexatic
phase in case of the clean system), and finally at a high temperature (deep
into the liquid state). The following points are worth mentioning:
The results for clean system in Fig. S1(a-c) are consistent with the KTHNY
picture of 2D-melting in a pure systems. Note that the particles in a hexatic
phase (Fig. S1(b)) preferentially move along the three principal directions of
the underlying triangular lattice of the solid.
The lowest temperature phase for the RP-system in Fig. S1(d) represents a
hexatic glass (or, an amorphous solid) in which the pinned impurities are
marked red, with no apparent positional order. The snapshot at
$\Gamma^{-1}=0.0140$ for RP-system in Fig.S1(e) appears inhomogeneous, where
the mobile particles close to pinned particles delocalize more. Other
particles tend to carry the motional signature of a hexatic phase in a clean
system. The snapshot in Fig.S1(f) for RP-system for $\Gamma^{-1}=0.0210$
indicates that delocalization of mobile particles is nearly engulfed the whole
space, except for small region surrounding the repulsive pinning centers.
The snapshots in the CP-system in Fig.S1(g-i) supports the notion of enhanced
solidity in the following manner: The solidity is nearly perfect at low $T$
(solidity here is stronger than in clean system) – the impurities are pinned
at commensurate locations of the underlying triangular lattice hold onto
perfect solidity. Unlike for RP-system, here the regions rich in pinning
centers anchors solidity around them, whereas, the regions relatively free of
impurity feature incipient melting. This is seen from Fig.S1(h) for
$\Gamma^{-1}=0.0182$, i.e. at a temperature just below the onset of melting.
As a result, the CP-system displays a local pockets of ‘remnant solidity’
around impurities. This remnant solidity weakens with $\Gamma^{-1}$, however,
persists up to a large temperature (Fig.S1(i) at $\Gamma^{-1}=0.0210$). Note
that the pinning centers maintain solid-like correlation at all $T$. The
remnant solidity has discernible effects on the bond-orientational
correlations $g_{\rm 6}(r)$, discussed in connection with Fig. 5(d) in the
main manuscript, and is elaborated further in the later part of this SM.

##  III Identification of critical distances of defects unbinding

Figure S2:  The distribution of distances between dislocation pairs: $P(r_{\rm
dp})$ shown for different phases. Panels (a-c), (d-f), and (g-i) show results
for clean, RP-, and CP-systems, respectively. $P(r_{\rm dp})$ at low-$T$ in
panels (a), (d), and (g) are sharply peaked at the typical distance of
dislocation pairs, though there are differences in details. The distribution
has a short tail for the clean system, and a very long tail for RP-systems,
and essentially no tail for CP-systems, indicating the enhanced solidity due
to commensurate pinning. The distributions for the three systems develop
smooth and healthy tails at a $T$, where defects unbinding sets in, as shown
in panels (b, e, h). Panels (c, f, i) show that they become indistinguishable
at large $T$, beyond the melting to a liquid.

In fig. S2, we report the distribution of the separation of the two
dislocations with equal and opposite Burger’s vector (here the two constitute
a bound pair of dislocations). We denote this distribution by $P(r_{\rm dp})$,
where $r_{\rm dp}$ is the separation in question. At lower $T$
($\Gamma^{-1}=0.0112$), Fig. S2(a) show that $P(r_{\rm dp})$ in a clean system
is sharply peaked at $r_{\rm dp}\approx a_{0}$, the average inter-particle
distance, with some additional weight for larger $r_{\rm dp}$. In contrast, a
weak yet very long tail in $P(r_{\rm dp})$ (Fig. S2(b)) is found for the RP-
system at low-$T$. This demonstrates the presence of unbound dislocations even
at the lowest $T$ ($\Gamma^{-1}=0.0056$), prohibiting true solidity. In
contrast, the first peak in $P(r_{\rm dp})$ for CP-systems at low-$T$ in Fig.
S2(g) is even sharper than the one in clean system (Fig. S2(a)). This reflects
an enhanced solidity in the CP-system. The nature of the $P(r_{\rm dp})$
changes with $T$, and when the dislocation pairs start unbinding, $P(r_{\rm
dp})$ begins to develops smooth tails, see Fig. S2(b, e, h). At these $T$’s
and beyond $P(r_{\rm dp})$ becomes qualitatively similar for clean, RP- and
CP-systems. The distributions become nearly identical deep into the liquid
phase (hight-$T$) as shown in Fig. S2(c, f, i).

Figure S3: Difference in $P(r_{\rm dp})$: The two panels display $\Delta
P(r_{\rm dp})$ in the pure system, taken between two pair of temperatures –
one in the solid and the other in the hexatic phase, both being close to
${}^{C}\Gamma^{-1}_{\rm SH}$. The corresponding traces of $P(r_{\rm dp})$ is
shown as the inset. The largest $r_{\rm dp}$ corresponding to a zero-crossing
is taken as $r_{dp}^{c}$.

Fig. S3 supplements the Fig. 3(d) in the main text, and illustrates the
insensitivity of $r_{dp}^{c}$ on the chosen $T$-values for its extraction, as
long as those temperatures are close to ${}^{C}\Gamma^{-1}_{\rm SH}$. We plot
in Fig. S3 the difference $\Delta P(r_{\rm dp})$ of $P(r_{\rm dp})$ for two
different pair of temperatures around the ${}^{C}\Gamma^{-1}_{\rm SH}$, while
the corresponding distributions are shown as the inset. Both plots yield
$r_{dp}^{c}=2.15a_{0}$ – the same value obtained in the main test for a yet
another pair of $T$’s.

Figure S4: integrated distribution (ID): as a function of $r_{\rm dp}$ for
different $Ts$. The attainment of unity of ID occurs sharply in the solid
phase and more gently beyond ${}^{C}\Gamma^{-1}_{\rm SH}$. The curvature of ID
is studied further to estimate the critical distance ($r_{dp}^{c}$), which is
presented as an inset.

An alternative extraction of $r_{\rm dp}^{c}$ for the pure system, employing
differently the same concept as above is the following: In Fig. S4, we plot
the integrated distribution (ID) of $P(r_{\rm dp})$ for various $T$’s. The
traces of ID shows a threshold behaviour marking the solid-to-hexatic
transition. ID comes up very steeply for the lower $T$’s corresponding the
solid phase (as inerred from the Fig. 2 of the main text), whereas the rise of
ID becomes distinctively gradual for the two higher $T$’s representing a
hexatic phase. In the inset, we plot the distribution of the change of the
slope of ID-traces in the solid phase, which attains the peak-value for
$r_{\rm dp}^{c}\approx 2.15a_{0}$.

Figure S5:  distribution of distances between disclination pairs: Panels (a-c)
represent the distribution of the distances of the disclination pairs
$P({\tilde{r}}_{\rm dp})$ at different $\Gamma^{-1}$. Panel (a) is the
$P({\tilde{r}}_{\rm dp})$ corresponds to the hexatic phase, which attains a
sharp peak around the lattice spacing reflecting the only presence of tightly
bounded disclinations in the phase. When the system hits the transition $T$,
the defects unbounds. The unbinding of these defects modifies the
distribution, as shown in panels (b) in the liquid phase, and remains
identical for the higher $\Gamma^{-1}$ value (panel (c)).

We also studied the thermal evolution of the distribution of the disclination
pair (${\tilde{r}}_{\rm dp}$) and we present our results in fig .S5 for three
representative $T$’s. This distribution $P({\tilde{r}}_{\rm dp})$ features a
single peak at ${\tilde{r}}_{\rm dp}\sim a_{0}$ for up to
${}^{C}\Gamma^{-1}_{\rm HL}$, and identify a second peak for larger $T$’s at
${\tilde{r}}_{\rm dp}\approx 2a_{0}$. Thus, we conclude that
${\tilde{r}}^{c}_{\rm dp}\sim r^{c}_{\rm dp}$.

##  IV Correlations:

The pair correlation is given by

$g(r)=\frac{1}{2\pi rN}\sum_{i=1}^{N}\sum_{j\neq
i=1}^{N}\langle\delta(r-\lvert\bf{r}_{i}-\bf{r}_{j}\rvert)\rangle,$ (1)

with $r=\lvert{\bf r}_{i}-{\bf r}_{j}\rvert$. We present the low-$T$ behavior
of $g(r)$ in Fig. S6(a). Its evolution in the clean (blue) and CP (green)
systems is nearly identical, consistent with the findings of Fig. 2(a) and
Fig. 2(b) in the main text. The $g(r)$ for the RP-system exhibit damped
modulations which prohibits solidity even at low temperatures. However, these
damped modulations are long-ranged than what’s expected in a liquid. Fig.
S6(b) shows $g(r)$ for large $T$ where all three systems turns into ‘liquid.’
A role reversal occurs – $g(r)$ for clean and RP-system overlap this time,
yielding a liquid-like behavior, while the long-range modulations for CP-
systems survives implying ‘remnant solidity’ as discussed already in other
context.

Figure S6: Pair correlation function: Panel (a) represents the pair
correlation (red: clean-, blue: random- and green: CP-system) for two
different $\Gamma^{-1}$ values. At the low $\Gamma^{-1}$ value, the clean
system’s behavior is comparable to the CP-system, while in the higher
$\Gamma^{-1}$, as shown in panel-b (for visual clarity, the traces of RP-
system is shown in a dotted line), it is comparable to the RP-system.  Figure
S7:  Orientational correlation function: of a clean system at different $Ts$.
The pure system displays the thermal evolution of $g_{6}(r)$ consistent with
KTHNY theory.

The $T$-dependence of orientational correlations in the clean system is shown
in Fig. S7 implies that $g_{6}(r)$ remains nearly independent of $\Gamma^{-1}$
for up to $\Gamma^{-1}=0.0145$. It shows a power-law decay $\sim
r^{-\eta(T)}$, in the hexatic phase for $0.0151\geq\Gamma^{-1}\geq 0.0162$
(exponents $\eta(T)$ are listed in Fig. S7). Finally, $g_{6}(r)$ features an
exponential decay in the liquid phase for $\Gamma^{-1}\geq 0.0165$. We find,
e.g., $g_{6}(r,\Gamma^{-1}=0.0140)\sim(r/a_{0})^{-0.003}$,
$g_{6}(r,\Gamma^{-1}=0.0.0157)\sim(r/a_{0})^{-0.01}$,
$g_{6}(r,\Gamma^{-1}=0.0159)\sim(r/a_{0})^{-0.23}$, while
$g_{6}(r,\Gamma^{-1}=0.0182)\sim{\rm exp}^{-0.24(r/a_{0})}$ and
$g_{6}(r,\Gamma^{-1}=0.0210)\sim{\rm exp}^{-0.49(r/a_{0})}$.

Figure S8:  Orientational correlation function, The $g_{\rm 6}(r)$ is
calculated in two separate components. (i) The dotted line shows the
contribution of $g_{\rm 6}(r)$ only for the particles which lie within a
distance of $2a_{0}$ from pinned particles at the higher $Ts$ beyond the
melting temperature and (ii) the solid traces shows the $g_{\rm 6}(r)$
contributed by all other particles which are more than $2a_{0}$ distance.

In order to address ‘remnant solidity’ in the CP-system at high $T$, we plot
the two separate components of $g_{\rm 6}(r)$ in Fig. S8. Once component
estimates the contribution of $g_{\rm 6}(r)$ (dotted trace) only for the
particles which lie within a distance of $2a_{0}$ from pinned particles. The
other component measures $g_{\rm 6}(r)$ (solid trace) contributed by all other
particles. We present data of such two contributions for two values of $T$
(both in the high-$T$ liquid phase). The message from Fig. S8 is clear – the
particles close to the pinning centers at commensurate locations of a perfect
crystal anchor the crystallinity around them, which is reflected in the
corresponding $g_{\rm 6}(r)$ which remains nearly constant at large $r$ like
in a solid. The particles further away from the pinned particles, on the other
hand, expectedly behaves like a liquid. Thus, it is the particles near pinning
centers which contribute to the ‘remnant solidity.’ As such, that pinned
particles must show a perfect solid like correlation by construction!


# missing

Entanglement Dynamics in
Quantum Continuous–Variable States

DOCTORAL THESIS

Ankit Kumar

Indian Institute of Technology Roorkee, India

_jointly supervised by_

Prof. P. Arumugam

Department of Physics

Indian Institute of Technology Roorkee

Roorkee 247667, India

Prof. Tomasz Paterek

Institute of Theoretical Physics and Astrophysics

Faculty of Mathematics, Physics and Informatics

University of Gdańsk, Gdańsk 80-308, Poland

July, 2023

©INDIAN INSTITUTE OF TECHNOLOGY ROORKEE, 2023
ALL RIGHTS RESERVED

headsepline=true Abstract

Due to the weakness of gravitational coupling, all quantum experiments up to
date in which gravity plays a role utilized the field of the Earth. Since this
field undergoes practically undetectable back-action from quantum particles,
it effectively admits a classical description as a fixed background Newtonian
field or spacetime. This argument strongly motivates theoretical and
experimental research towards a demonstration of gravitation between two
quantum masses, as this is one of the most straightforward scenarios where
quantum features of gravity could be observed. Several proposals studied the
possibility of generating entanglement between two massive objects. Along the
same lines, with a particular focus on gravity, this thesis introduces general
tools to tackle interaction-mediated entanglement and applies them to two
particles prepared in continuous-variable states.

In order to pursue this aim systematically, this dissertation begins by
introducing methods to precisely simulate the dynamics of quantum systems
coupled by weak interactions. We improve the accuracy of the numerical
implementation of Cayley’s operator and develop a methodology to avoid
reflections from numerical infinities. We derive a condition under which a
product state from the laboratory (LAB) perspective remains a product state in
the center-of-mass (COM) frame, which reveals that only certain states are
transformed into disentangled states. Even though the primary focus is on
gravity, all the developed methods apply to arbitrary central interactions,
and considerable parts of this thesis are devoted to explicit demonstrations
of this versatility. Accordingly, the first application is to investigate the
head-on collision in the Rutherford experiment, with the projectile treated as
(realistic) localized wave packets. We observe various nonclassical effects in
the average trajectories and trace them back to the convexity properties of
the Coulomb potential with the help of Jensen’s inequality. The concluding
chapter also comments on the projectile-target entanglement.

Our next goal is to simplify the possible observation of weak gravitational
entanglement in an inevitably noisy laboratory. The basic idea is to amplify
correlations by pushing the particles toward each other, hoping that an ever-
increasing gravitational interaction will automatically lead to a higher
accumulated entanglement. A toolbox is developed that quantifies the
entanglement gain between the two particles directly in the COM frame of
reference, thereby eliminating the need for inverse transformations back to
the LAB frame. We start with the standard practice of the second-order
truncation of quantum Newtonian potential, which has long kept the
mathematical complexities at a minimum by forcefully constraining the system
into the regime of (very well-understood) Gaussian Quantum Information. While
it is known that an analytical solution exists, we utilize Ehrenfest’s theorem
to derive the covariance matrix in an exact closed form. The resultant
entanglement is insensitive to relative motion between the two particles.

The less-understood non-Gaussian regime triggered by the cubic and higher-
order potentials is considered next. We develop a hybrid analytical-numerical
scheme to faithfully estimate the entanglement gain with the help of
algorithms in Google TensorNetwork. The entanglement is found to be sensitive
to relative motion only when the system evolves into the non-Gaussian regime.
We prove that the position-momentum correlations originate from the force
gradient in relative motion. A derivation of closed forms for the non-Gaussian
entanglement gain follows with informed guesswork. In experiments, it will be
challenging to screen the system from all interactions but gravity. With this
in mind, we develop tools to quantify the entanglement with multiple central
forces acting simultaneously.

As the final application, the thesis discusses an entanglement-based test of
the Modified Newtonian Dynamics (MOND), a candidate explanation of dark matter
effects which proposes to modify Newton’s second law and/or the gravitational
force law for accelerations smaller than $\sim 10^{-10}$ m/s2. One verifies
that the masses recently cooled by the Aspelmeyer group in Vienna, when
separated by a distance of a few times their radius, are into the regime of
accelerations where MOND is relevant. Accordingly, the tools developed in this
thesis offer an opportunity to test the assumptions behind MOND through
entanglement between two nearby quantum masses. We develop an experiment where
departures from Newtonian gravity are certified by simply witnessing the
entanglement generation starting from thermal states. Publications

## Contributions in Peer-Reviewed Journals

  1. 1)

Nonclassical Trajectories in Head-On Collisions
Ankit Kumar, Tanjung Krisnanda, P. Arumugam, and Tomasz Paterek
Quantum 5, 506 (2021) arXiv:2011.06470

  2. 2)

Continuous-Variable Entanglement through Central Forces: Application to
Gravity between Quantum Masses
Ankit Kumar, Tanjung Krisnanda, P. Arumugam, and Tomasz Paterek
Quantum 7, 1008 (2023) arXiv:2206.12897

  3. 3)

Probing Modified Gravity with Entanglement of Microspheres
Ankit Kumar, Yen-Kheng Lim, P. Arumugam, Tom Złośnik, and Tomasz Paterek
Physical Review D 109, L101501 (2024) arXiv:2306.14938

## Contributions in Peer-Reviewed Conference Proceedings

  1. 1)

Closest Approach of a Quantum Projectile
Ankit Kumar, Tanjung Krisnanda, P. Arumugam, and Tomasz Paterek
Journal of Physics: Conference Series 1850, 012074 (2021) arXiv:2112.13296

## Contributions in Public Scientific Repositories

  1. 1)

An Accurate Pentadiagonal Matrix Solution for the Time-Dependent Schrödinger
Equation
Ankit Kumar, and P. Arumugam
Zenodo.7275667 GitHub/vyason/Cayley-TDSE arXiv:2205.13467

###### Contents

  1. Contributions in Public Scientific Repositories