# Non-Normal Interactions Create Socio-Economic Bubbles

Didier Sornette  Sandro Claudio Lera  Jianhong Lin  Ke Wu  Institute of Risk
Analysis, Prediction and Management (Risks-X), Academy for Advanced
Interdisciplinary Sciences, Southern University of Science and Technology,
Shenzhen, China  Connection Science, Massachusetts Institute of Technology,
Cambridge, USA Blockchain and Distributed Ledger Technologies, Institute of
Informatics, University of Zürich, Andreasstrasse 15, 8050 Zürich, Switzerland
UZH Blockchain Center, University of Zürich, Andreasstrasse 15, 8050 Zürich,
Switzerland

###### Abstract

In social networks, bursts of activity often result from the imitative
behavior between interacting agents. The Ising model, along with its variants
in the social sciences, serves as a foundational framework to explain these
phenomena through its critical properties. We propose an alternative generic
mechanism for the emergence of collective exuberance within a broad class of
agent-based models. We show that our model does not require the fine-tuning to
a critical point, as is commonly done to explain bursts of activity using the
Ising model and its variants. Instead, our approach hinges on the intrinsic
non-symmetric and hierarchical organization of socio-economic networks. These
non-normal networks exhibit transient and unsustainable surges in herd
behavior across a wide range of control parameters even in the subcritical
regime, thereby eliminating the need for the - arguably artificial - fine-
tuning proximity to a critical point. To empirically validate our framework,
we examine the behavior of meme stocks and establish a direct linkage between
the size of financial bubbles and the degree of non-normality in the network,
as quantified by the Kreiss constant. Our proposed mechanism presents an
alternative that is more general than prevailing conceptions of instabilities
in diverse social systems.

††journal:

## Introduction

Many complex dynamical systems are characterized by periods of relative
stability which are interrupted by transient regimes during which the dynamics
exhibits sudden bursts of unsustainable growth (a “bubble”) or shifts suddenly
to another attractor. A large corpus of knowledge and methods have been
developed in the last two decades to account for these phenomena, which are
based on the underlying concept of tipping points, wherein a critical
threshold is reached at which the system bifurcates to a new state [1, 2, 3,
4, 5]. Typical financial models that account for such instabilities contain
two classes of traders: fundamentalists who maximize their expected utility
function and noise traders [6, 7]. Noise traders are usually assumed to
influence each other according to an Ising-like dynamics, with interaction
dependencies captured by an adjacency matrix $\mathbf{A}$ and interaction
strength captured by a coupling constant $\kappa$. When the imitation strength
between noise traders is large enough, collective social behavior can occur,
such as polarization of noise traders toward buying (selling), which in turn
creates bubbles (crashes) [8, 9].

These bubbles and crashes are generally associated with the underlying Ising
phase transition separating a disordered opinion regime for low imitation
strength $\kappa$ from an ordered regime where all noise traders tend to be
synchronized. In all existing models of this type, bubbles requires the
imitation strength $\kappa$ to be close to or larger than a critical value
$\kappa_{c}$ associated with the underlying phase transition. In other words,
in this class of models, bubbles and crashes are the signatures of the fact
that the financial market has entered a “critical regime”, in the technical
sense of the emergence of collective order in the decisions of a large
fraction of traders. There is an extensive literature on agent-based models
and generalized Ising models [10, 11, 12, 13, 14]. To the best of our
knowledge - in all cases - the abnormal stylized facts, such as excess
volatility and transient bubbles and crashes, require the system to be close
to, at or slightly above the critical point in the ordered polarizing regime.

In this work, we document a general mechanism for the nucleation and growth of
transient bubbles. We suggest that this new mechanism is much more general
than previously existing ones and likely to be often the dominant process at
work, because it does not require the fine-tuning close to or sweeping [15] of
the system over a critical point. Our proposed mechanism is based on the fact
that social influence is characterized by two critical properties: it is
directed and hierarchical [16, 17]. For instance, in our example of financial
markets, the influence of a famous investor on a retail investor is likely
much larger than the other way around. Together, these two ingredients give
rise to networks with non-normal adjacency matrices $\mathbf{A}$, called non-
normal networks [18, 19, 20], whose associated dynamical systems are known to
induce transient bursts [21, 22, 23, 24]. Interpreted in terms of socio-
economic interactions, these transient bursts are responsible for short-lived
social contagion even well below any critical threshold. We demonstrate this
mechanism in the context of the formation of financial bubbles and their
following crashes that lead to enormous economic losses. Analyzing Reddit
discussion forums of meme stocks, we show that patterns of influence are
highly non-normal, and that the rate of reciprocity is dependent on a user
hierarchy. Using a previously validated agent-based model [25, 26], we show
that non-normal networks give rise to transient bubbles even when the
imitation strength is sub-critical. Intuitively, some traders are more
influential than others and information does not spread evenly but along
cascading circuits. Our work thus provides a qualitative proposal that
financial systems are intrinsically generating crises [27]. Due to the broad
applicability of stochastic interaction models [28], our model is expected to
generalize in a straightforward manner to other hierarchical socio-economic
systems, and help explain wide-spread phenomena such as social bubbles [29]
and herding in opinion dynamics [11].

## Results and Discussion

### Agent Based Price Simulation

Agent-based models (ABMs) have become a popular tool in interdisciplinary
research over the last decades [30, 31, 32], primarily due to their
flexibility in accounting for heterogeneous and non-linear interactions. Here,
we implement an ABM that simulates a financial market consisting of
fundamentalists and noise traders who trade a risky and a risk-free asset [25,
26]. The risky asset is a dividend paying stock. The risk-free asset pays a
constant return in each time-step and represents a bank account or risk-free
government bond. Each trader formulates their excess demand for the next time
step and the price $P_{t}$ of the risky asset is calculated as the Walras
equilibrium in which supply equals demand. Fundamentalists are rational risk-
averse investors who invest by maximizing their expected utility under a
constant relative risk aversion utility function. At each time step, they
allocate their wealth between the two assets in order to maximize their
expected utility over the next period. They buy (sell) or sell (buy) the risky
(riskless) asset in order to move their portfolio towards their desired
optimal portfolio, given the information they have on price and volatility of
the risky asset. They play a role akin to price-setters, in contrast with
noise traders who are pure price-takers. Their demand for the risky asset
combines with that of the noise traders to fix the price at each time step. We
refer to the Supplementary Notes for more details on the fundamentalists. 111
The Supplementary Notes are available from the authors upon request.

The crucial component for our model are the noise traders. Their investment
strategy is based on an Ising-like social influence model, where they can be
modeled as nodes in a network with directed edges. While previous
implementations [25, 26] have considered a mean-field approximation in terms
of a fully connected, symmetric network, here we make the network topology an
explicit component of our model. As will become clear below, relaxing the
assumption of fully connected, symmetric interactions is not only realistic,
but also a crucial ingredient for the observation of sub-critical bubbles.
Each trader $i$ is considered to be in one of the two possible states, $+1$
(the noise trader holds the risky asset), and $-1$ (the noise trader holds the
risk-free asset). The states are denoted as $s^{i}=\pm{1}$, respectively. The
directed network of $N$ nodes (noise traders) is described by its adjacency
matrix $\mathbf{A}=\\{a_{ij}\\}$, where $a_{ij}=1$ if there exists a directed
edge (influence) from node $j$ to node $i$, and $a_{ij}=0$ otherwise.
Introducing further the transition matrix $\mathbf{M}$ which is proportional
to $\mathbf{A}$ (see Methods section for details) allows us to represent the
state dynamics of the Ising-like noise traders as

$\triangle{\vec{s}(t)}=\mathbf{M}\vec{s}(t)$ (1)

At each time step $t$, the collective opinion (“magnetization” in the Ising
language) $m_{t}$ is defined as

$m_{t}=\frac{1}{N}\sum_{i=1}^{N}s^{i}_{t}~{}~{}~{}\in[-1,1].$ (2)

This system remains stable as long as the real parts of all eigenvalues of
$\mathbf{M}$ are negative. As is well known, by continuously tuning
$\mathbf{M}$, systems whose linear stability is controlled by (1) can undergo
a bifurcation (or phase transition) from a stable fixed point with zero
average change of spin to a state where all spin change state to align to each
other (a state described by higher-order terms beyond the linear expansion
$\mathbf{M}\vec{s}(t)$). The existence of such states has been related to the
emergence of financial bubbles (crashes), diagnosed by the existence of
transient super-exponential growth (loss) [33, 9]. System (1) is well-known to
be a generic representation of the aggregate behavioral outcomes of financial
markets where investors are driven by group psychology and herding [6, 34, 35,
36]. Moreover, as will become clear below, our results are not restricted to
the detailed assumptions of our model and generalize to a large class of
stochastic interaction models (from hereon referred to as ‘Ising-like’). In
the remainder of this article, we will analyze the sub-critical regime of
$\mathbf{M}$, but with $\mathbf{M}$ being non-normal. We will show that this
non-normal structure gives rise to transient dynamics that induce bubbles and
crashes much like at or close to criticality in normal networks.

Figure 1:  Difference between normal and non-normal transients. (a) and (b)
show the geometry of pseudo spectra. In each plot, the contours represent the
boundary of $\sigma_{\epsilon}(\mathbf{M})$ for two values of $\epsilon$ (see
Methods). (c) Initial, transient and asymptotic behavior of
$\left|\left|e^{t\mathbf{M}}\right|\right|$ for a non-normal matrix
$\mathbf{M}$. The graph indicates that here $\sigma(\mathbf{M})<0$ and hence
that the asymptotic behavior is stable. However, the asymptotic behavior is
not at all predictive of the transient behavior in case $\mathbf{M}$ is non-
normal. Plot (d) shows the eigenvalues (black dots) and some
$\epsilon$-pseudospectra for a normal matrix (different colors represent
different values of $\epsilon$). All eigenvalues as well as the epsilon-
spectral lines are confined to the left half plane of $\mathbb{C}$.
Accordingly, $||s_{t}||$ decays exponentially as shown in the inset plot. Plot
(e) shows the case of a non-normal matrix $\mathbf{M}$. While its eigenvalues
are also confined to the left half plane of $\mathbb{C}$, its $\epsilon$
spectral lines are not. According to inequality (3), in the inset plot, we see
intermittent transient growth before the asymptotic decay as in Plot (c).

### Non-Normal Matrices

Early contributions to the study of non-normal matrices have originated from
hydrodynamics, where non-normality plays a role in the emergence of turbulence
[21]. Ever since, non-normal matrix theory has helped explain phenomena such
as perturbations in ecosystems [37], non-Hermitian quantum mechanics [38],
population dynamics [39], synchronization of optoelectronic oscillators [40],
amplification of neural activities [23], chemical reactions [41], network
synchronization [18, 42, 43], neuronal networks [44, 45], ecological reactions
[46], and the formation of Turing patterns [24, 47].

A matrix $\mathbf{M}$ is called non-normal if
$\mathbf{M}^{T}\mathbf{M}\neq\mathbf{M}\mathbf{M}^{T}$. Non-normal matrices,
unlike their normal counterparts, contain at least one set of eigenvectors
which are non-orthogonal to each other. We refer to the Methods section for a
primer on non-normal matrices. For the purposes at hand, one key property of
non-normal matrices $\mathbf{M}$ is that the dynamics they describe via (1) is
significantly different from the long-term asymptotic behavior governed by the
largest eigenvalue. It is well-known that the asymptotic behavior for
$t\to\infty$ is governed by the largest real-part of all eigenvalues of
$\mathbf{M}$. In particular, so long as all real parts of $\mathbf{M}$ are
negative, the dynamics is stable. Our main interest lies in transients, which
are described by intermediate values of $t$.

It can be shown that, for intermediate times, there is transient growth
according to

$\sup_{t\geqslant
0}\left|\left|e^{t\mathbf{M}}\right|\right|\geqslant\mathcal{K}(\mathbf{M}).$
(3)

Importantly for our application below, this initial transient growth is
exponential (Figure 1(c)). Given an interaction matrix $\mathbf{M}$ as in (1),
we can calculate the Kreiss constant $\mathcal{K}(\mathbf{M})$ (see Methods)
to obtain lower bounds for the transient growth of net magnetization. An
example of such transient growth is shown in Figure 1(c) and in the inset of
Plot (e). As we shall see, these transients are responsible for socio-economic
bubbles in our agent-based model.

### Level-Dependent Reciprocal Connections

A system such as eq. (1) with non-normal $\mathbf{M}$ can be interpreted as a
dynamical process on a complex asymmetric and partly hierarchical network.
Such non-normal networks have been observed in a wide variety of biological
and socio-economic networks [48, 19, 20], and their role in the transmission
of noise has been studied [49]. Recall that asymmetry of $\mathbf{M}$ is a
necessary, but not a sufficient condition for non-normality. For instance,
consider a simple cyclical, directed network of three nodes $\\{X,Y,Z\\}$
where $X\to Y,Y\to Z$ and $Z\to X$. An adjacency matrix with such cyclical
symmetry still gives rise to a normal adjacency matrix. The condition
$\mathbf{M}\mathbf{M}^{T}\neq\mathbf{M}^{T}\mathbf{M}$ is instead satisfied
when the network is directed and hierarchical, which are both intrinsic
properties of socio-economic systems [17, 50, 51]. Indeed, it is well-
established in anthropology that people organize their social interactions
along one or more of the four elementary modes: communal sharing, authority
ranking, equality matching, and market pricing, with the latest one being
particularly relevant in finance [52, 53]. Of these four modes, only equity
matching is symmetrically directed, while the others exhibit some direction
asymmetry and a degree of hierarchy. Assuming that non-normal matrices govern
the interaction of humans in socio-economic and financial systems is thus
natural.

Recently, methods to generate non-normal networks have been proposed by taking
into consideration asymmetrical reciprocity [48] and hierarchy [19] that are
typical of non-normal systems. Drawing from these insights, we implement here
an algorithm that allows us to control the rate of non-normality along with
the number of top nodes, that can be interpreted as thought leaders. In
contrast to previous work [48], our rate of reciprocity explicitly depends on
the hierarchical level which is a realistic addition as has been evidenced in
a variety of biological and social non-normal networks [19].

The non-normal network with a total of $N$ nodes is initialized with $N_{0}$
so-called top nodes. These top nodes account for the largely independent
$N_{0}$ backbones of the communication network common to typical
hierarchically, non-normal networks [48, 19]. The remaining $N-N_{0}$ nodes
are added to the existing network sequentially, one node at a time. Each newly
added node receives $m$ in-edges, i.e. channels of communication through which
it can be influenced. The source of each such edge is selected with
probability proportional to the existing nodes’ out-degree. As is well-known,
this type of preferential attachment creates a skewed degree distribution
whereby the network is dominated by a few central nodes [54, 55]. Once the $m$
source nodes of an added node are determined, each of the $m$ newly formed
directed edges may be reciprocated with some probability $\rho$. Figure 2(a-c)
exemplifies several cases where $\rho$ is equal to a fixed value $\theta$.
Setting $\theta\ll 1$ gives rise to strongly non-normal systems [48], whereas
$\theta=1$ recovers a symmetric, i.e. normal, network.

Figure 2:  Hierarchical levels and rates of reciprocity in complex networks.
(a-c) Three examples of directed networks with different levels of non-
normality, different number of top nodes, and different but constant rates of
reciprocity $\rho$ given by different values of $\theta$. The colors indicate
the hierarchical (in general non-integer) level $\ell$ of the nodes. A value
of $\theta=0$ means no edge can be reciprocated, whereas a value of $\theta=1$
means every edge is reciprocated (since level-dependence is ignored). Network
(c) is then normal, since it is symmetric. (d) Empirical analysis of the
Blackberry subreddit network. The rate of reciprocity is not constant, but a
function of the hierarchical level, $\rho=\rho(\ell)$. The higher the level,
the higher the rate of reciprocity, up to some level of saturation. Error bars
represent standard deviations of average rate of reciprocity measured across
all users. The red line shows the sigmoid function that best fits the data.

Based on empirical evidence [19] and to reflect the fact that nodes that are
higher up in the hierarchy are harder to be influenced, we assume further that
this probability $\rho$ is modulated by the hierarchical level $\ell$ of each
node, $\rho=\rho(\ell)$. The lower the node in the hierarchy (the larger
$\ell$), the more likely the node is reciprocated. Loosely speaking, the
hierarchical level $\ell$ of any node $i$ is defined as the shortest path from
a top-node to node $i$. More precisely, the level $\ell$ is defined as the
trophic hierarchical level [56, 57]. We have analyzed the Reddit discussion
forum for the Blackberry meme stock below. An edge is drawn from user $i$ to
user $j$ if $j$ replies to a comment of user $i$. In Figure 2(d), we show that
the rate of reciprocity is not constant, but an increasing function of $\ell$.
In other words, the more popular a user’s comments, the less likely that user
is to reciprocate (comment on) any given edge. For social communication
networks, this observation is natural and rationalized as the approximately
constant finite capacity of any given individual to respond to comments.

We detail in the Methods section our algorithm to generate non-normal
adjacency matrices. Such a matrix $\mathbf{A}$ has six parameters:
$N,N_{0},m,\theta,a$ and $b$. The four parameters $m,N,a$ and $b$ play a
subordinate role in the qualitative interpretation of our results. For the
remainder of this paper, we thus fix $N=1000$, $m=2$, $a=2.552$ and $b=3.668$
unless mentioned otherwise. The parameters $\theta$ and $N_{0}$, on the other
hand, have qualitatively important implications on the behavior of our model.
The parameter $\theta$ characterizes the asymmetric nature of the system. The
smaller $\theta$, the more directed the network, and the less the top nodes
may be influenced. The parameter $N_{0}$ denotes the number of top nodes. If
$\theta$ is small, then $N_{0}$ may be interpreted as the number of (largely)
independent, leading opinions in the system. As we will now show, our results
do not depend on the specific implementation of the non-normal mechanism, but
are primarily mediated by the Kreiss constant (see Primer on Non-Normality in
the Methods section).

### Transient Bubbles Induced by Non-Normal Interactions

Figure 3:  Normal vs. non-normal price dynamics. (a) Magnetization (2) for an
ABM simulation with a non-normal interaction matrix $\mathbf{M}$ for different
parameter constellations. The net magnetization of the $N_{0}$ opinion leaders
is shown in grey lines in the background. The opinion leaders change their
opinion randomly at rate $p^{\pm}/2$, as dictated by equation (4). (b) Same as
in (a), but for a symmetrized (i.e. normal) interaction matrix
$\frac{1}{2}\left(\mathbf{M}+\mathbf{M}^{T}\right)$. The parameter $\kappa$ is
chosen at a sub-critical value, hence the net magnetization is, on average,
$0$. In contrast to (a), the transient bursts are much less pronounced. (c)
Price trajectory generated by the agent-based model dynamics with
magnetization from (a) and (b), respectively. Only the non-normal matrices
induce bubbles. (d) Bubble size as a function of Kreiss constant for different
parameter constellations. (e) Bubble steepness as a function of numerical
abscissa for different parameter constellations. (f) Bubble size as a function
of fraction of nodes receptive to antagonistic opinion. Error bars represent
standard deviations obtained from 100 simulations.

Building on the above insights, we now run agent-based simulations with a non-
normal adjacency matrix $\mathbf{M}$. While it has been well-established that
the formation of bubbles in agent-based models with Ising-like interactions is
associated with the proximity to the critical point of the underlying Ising
model with the occurrence non-zero net opinion (non-zero magnetization) [58,
59, 9], we investigate here the regime where the net magnetization (net
opinion) fluctuates around zero (sub-critical phase). In the following
analysis, we therefore set the coupling strength $\kappa$ to a sub-critical
value $(\kappa=0.98)$. The eight panels of Figure 3(a,b) confirms that $m_{t}$
fluctuates around zero approximately symmetrically, as expected from the fact
that the imitation strength $\kappa$ has been chosen so that the underlying
Ising model is subcritical. Furthermore, for fixed parameters
$(N_{0},\kappa)$, we compare two types of social networks: $\theta=0$ and
$\theta=1$. The former corresponds to a case of no reciprocity, and hence
large non-normality of $\mathbf{A}$ and hence $\mathbf{M}$. The later
corresponds to an almost symmetric - and hence normal - matrix $\mathbf{A}$,
which coincides with a much less non-normal matrix $\mathbf{M}$ (see also
Supplementary Notes). Comparing Figure 3 (a) and (b), we see that the strongly
non-normal case (small $\theta$) corresponds to much more pronounced long-
lived deviations of the magnetization from its zero average, as is expected
from transient dynamics (Figure 1(c) and (e)).

The right most column of Figure 3 shows the associated price $P_{t}$ as a
function of time, obtained as the Walras’ equilibrium between fundamentalists
and noise traders (Supplementary Note 1). It is striking to observe the
drastic differences in the price dynamics between highly non-reciprocal (non-
normal) interactions compared to reciprocated (normal) ones. In the former,
very strong price peaks are preceded by periods of strong price growth,
following by fast large asymmetric drawdowns. This qualifies the existence of
large amplitude bubbles as a clear diagnostic of this type of non-normal
networks. In contrast, for normal networks (and weakly non-normal $\mathbf{M}$
matrices), the price dynamics appears compatible with a standard geometric
Brownian motion at least at the qualitative level. As we now show, both of
these phenomena are explained as a function of the transients induced by non-
normality.

A hallmark of a financial bubble is the existence of unsustainable super-
exponential price growth [9, 33, 60]. The term super-exponential refers to a
special regime where the price grows much faster than exponential, with its
growth rate growing itself as a function of time, compared with exponential
growth characterized by a constant growth rate. Within our ABM, it can be
shown [25] that, to a first approximation, the price is an exponential
function of the net magnetization, $P_{t}=Ce^{cm_{t}},$ where the scaling
coefficient $c>0$ is a function of the model parameters. This result is a
direct consequence of the fact that, on short time-scales, the price is mostly
driven by the demand of noise-traders. Thus, our results are not specific to
the details of our implementation of the social imitation model and are
expected to hold for general classes of stochastic interaction systems [61]
and common assumptions about fundamentalists [35].

We recall that $m_{t}$ is defined as the average state across all trader
states $\vec{s}(t)$. The states $\vec{s}(t)$ are themselves governed by
equation (1) involving the non-normal interaction matrix $\mathbf{M}$, such
that $\vec{s}(t)\sim e^{\mathbf{M}t}\vec{s}(0)$. In a globally stable regime,
all eigenvalues of $\mathbf{M}$ associated with the stable equilibrium
$\vec{s}(t)=0$ are negative (Supplementary Note 1). The standard expectation
is thus that $m_{t}$ remains small and thus $P_{t}$ should not exhibit
abnormal fluctuations. But this is forgetting the transients induced by the
non-normality of $\mathbf{M}$. As shown in Figure 1(c,e), the asymptotically
stable fixed-point $\vec{s}=0$ is punctuated by initially exponential,
repelling dynamics over finite time scales. Furthermore, inequality (3)
provides us with a lower bound of the size of these transients, which are
mainly a function of the Kreiss constant $\mathcal{K}(\mathbf{M})$. Combining
$P_{t}=Ce^{cm_{t}}$ with the transient approximately exponential growth of
$m_{t}$, we thus predict the occurrence of finite lived bubbles qualified as
transient super-exponential growth of price. Here, the super-exponential
behavior is approximately described by an exponential of an exponential.

This reasoning leads us to hypothesize that various properties of emergent
bubbles should be mostly driven by metrics characterizing the non-normal
nature of the network of interacting noise traders. The most obvious candidate
to characterize non-normal networks being the Kreiss constant
$\mathcal{K}(\mathbf{M})$, we hypothesize that the dependence of the size of
emerging bubbles on parameters $(N_{0},\theta,\kappa)$ should reduce to a sole
dependence on $\mathcal{K}(\mathbf{M})$. To test this dependence, we measure
the size of the bubbles as the difference in price between the beginning and
the end of a regime of super-exponential growth (see Supplementary Note 5 for
details). It turns out that the sizes of bubbles are almost entirely explained
by just the Kreiss constant, for various parameter constellations, as shown in
Figure 3.

More precisely, for different parameter combinations of
$(N_{0},\theta,\kappa)$, we generate $100$ price simulations according to the
following procedure. We first generate a matrix $\mathbf{M}$ according to our
algorithm described above, and we subsequently run the agent-based model to
generate a time-series with $25,000$ time-steps, corresponding to $100$ years
(considering 250 trading days per calendar year). On each time series, we
measure the size of all bubbles. These sizes are subsequently averaged across
all 100 simulations, with the standard deviation serving as error bars. Figure
3(d) demonstrates the existence of a collapse of all curves when the average
bubble sizes are plotted as a function of the Kreiss constant
$\mathcal{K}(\mathbf{M})$ of the non-normal matrix $\mathbf{M}$. A large
Kreiss constant is associated with large bubble sizes, for different network
non-normality and social coupling $\kappa$.

The theory of transients does not only make a prediction about the size of the
transients, but also about their steepness. As visualized in Figure 1(c), we
expect the steepness of the transients, and therefore of the magnetization and
then price, to be increasing in the numerical abscissa $\omega(\mathbf{M})$.
Figure 3(e) shows that the bubble steepness is indeed mostly explained by the
numerical abscissa $\omega(\mathbf{M})$ of the non-normal network, for
different sets of parameters.

These results underline the crucial role of the non-normality of the noise
traders interaction network in shaping the bubbles and their properties in our
agent-based model. The key insight is that different parametrization indeed
all collapse onto this scaling law relating the bubble size to the Kreiss
constant and the bubble steepness to the numerical abscissa of the matrix
$\mathbf{M}$ characterising the social network of noise traders. The intuition
behind this mechanism is related to the hierarchical organization associated
with non-normal matrices [19]. The information flow is structured such that it
is driven by a few leaders, thereby reducing the effective number of
independent opinions and increasing the level of polarization.

Finally, we test the effect of an influential contrary opinion on the size of
the bubble. We first simulate a price dynamics with a single top node,
$N_{0}=1$, with $N=1000$, $m=2$ and $\theta=0$. Upon formation of a bubble
(continuous price increase for 50 time-steps), we interfere with the system by
holding a contrary opinion (opposite state of $N_{0}$) for $\Delta t$ time-
steps. That contrary node is connected to a fraction $f$ of all nodes with the
exception of the top node. Figure 3(f) shows the dependence of the bubble size
as a function of $f$. The larger $f$ is, the larger is the decrease of the
bubble size. This result is encouraging, suggesting that an external
controller can partially counteract the development of sub-critical bubbles.
However, this control remains limited even when $f\to 1$ and requires one to
reach to a large fraction of noise traders. A more scalable approach, that we
leave for future research, would be a minimal influence of a few key nodes in
order to achieve overall noise-cancellation, as has been recently shown in
communication networks [49], and which could improve on the more standard
market intervention involving large balance sheet build-up of major financial
agents such as a central bank [62].

### Non-Normal Communication in Meme Stock Trading

Figure 4:  Evolution of real and simulated meme stock prices based on
discussions on Reddit forums forming non-normal networks of influence. (a) Co-
evolution of Reddit network Kreiss constant and Blackberry meme stock price.
(b) ABM price time-series resulting from the simulated magnetization $m_{t}$
of noise traders, with distinct peak around the same time as the real price
trajectory. The bands around the thick lines represent confidence intervals
obtained across 100 simulations.

So far, our assessment of bubbles has relied on agent-based simulations, where
we can control the experimental conditions. The difficulty with empirical data
is that, in general, one cannot observe the matrix $\mathbf{A}$ that governs
trader interactions. On social trading platforms, such as eToro, interactions
can be measured precisely [55], but the trading volume relative to the entire
market is small, such that its influence on the price is negligible. This is
not the case for so-called meme stocks which have enjoyed recent popularity.
Driven primarily by retail traders, meme stock trading activity has been shown
to be largely influenced by Reddit discussion forums [63, 64, 65].

Reddit is organized into subreddits on which specific topics are being
discussed. Users interact by submitting new posts and adding comments to
existing posts or comments. Here, we analyze the posts and comments related to
four popular meme stocks (Blackberry, Nokia, GameStop and AMC) under the
famous subreddit r/wallstreetbets (also known as WallStreetBets or WSB) that
has become notable for its colorful and profane jargon, aggressive trading
strategies, and for playing a major role in the GameStop short squeeze in
early 2021. For each of the four stock, at time $t$, we draw a directed edge
from user $J$ to user $K$ if $K$ has commented or replied a stock-related text
by user $J$ in the time interval $[t-\Delta t,t]$. In other words, $K$ has
been influenced by $J$’s action in the past $\Delta t$ days. With this
procedure, for each meme stock, we extract a dynamically evolving influence
network $\mathbf{A}(t)$, for which we can measure the Kreiss constant
$\mathcal{K}(t)$.

The evolution of the Kreiss constant, along with the trading price, is shown
in the top plot of Figure 4 for the Blackberry stock (see Supplementary Note 6
for similar plots for the other three meme stocks). The two most prominent
price peaks around January 2021 and June 2021 coincide with the two largest
peaks of the Kreiss constant trajectory. This gives force to our proposition
that increased non-normality (quantified by large values of the Kreiss
constant) favours the occurrence of transient super-exponential price behavior
(bubbles) associated with the transient growth of perturbations before their
relaxation. In other words, we interpret the presence of financial bubbles in
these meme stocks as reflecting at least partially the asymmetric hierarchical
structure of the reddit discussion forum that induced a polarized bullish
opinion among retail traders, which then in-turn pushed the price up. And
indeed, the mostly mentioned words in Jan 2021 among the reddit submissions
and comments related to BlackBerry are rocket, bb, gme, shares, and buy. In
particular, among the 58,793 mentions of “rocket” and 12,132 mentions of “buy”
in 2021, 69% of “rocket” and 43% of “buy” were in January. More generally, as
is shown in Supplementary Note 6, we do find a positive correlation between
the Kreiss constant and price bubbles across meme stocks. This suggests that
the non-normal structure of the Reddit meme stock discussion forum is an
important driver of the observed price instabilities.

A strong asymmetric hierarchical structure of the reddit discussion forum
quantified by a large value of the Kreiss constant provides a powerful
catalysis for the emergence of transient price bubbles. But in a world of mass
communication and a plurality of social media and information channels, not
all observed price perturbations can realistically be expected to be
attributed to the discussions on Reddit forums. As for general dynamics with
non-normal operators, not all perturbations go through a non-monotonous
transient amplification. The realized trajectory of the transient very much
depends on the projection of the perturbations onto the pseudo-eigenvectors
[22, 24, 66]. Not all large Kreiss constant values should thus lead to a
bubble, as the market dynamics is more complex and cannot just be reduced to
one source of influence. The mapping between large Kreiss constant and bubbles
becomes rigorous when considered in terms of ensembles of price trajectories.
To show this, and to isolate the effect of discussions on Reddit, we insert
the empirical Blackberry discussion forum influence network $\mathbf{A}(t)$ as
input to our ABM. We simulate the resulting price time-series $100$ times and
keep track of the average price as well as its average standard deviation (see
Supplementary Note 6 for details). In the bottom plot of Figure 4, one can
observe indeed that price spikes coincide - in their ensemble average - with
peaks of the Kreiss constant, supporting our hypothesis that the non-normality
in the Reddit discussion forum contributes to explain the observed price
bubbles. While meme stocks happen to be a prime example of price dynamics that
are driven by opinion leaders [64], human interactions in general and
financial markets in particular are known to be subject to non-normal
interactions [52, 50, 36]. We thus hypothesize that the mechanism described in
this article is an important contribution to the manifold observations of
bubbles in socio-economic systems [58, 29, 33].

Our simulations exemplify the possibility to diagnose regimes of financial
instabilities by measuring the evolution of the Kreiss constant of the
underlying network of social interactions between traders. Periods in which
the Kreiss constant is large should be interpreted as prone to bubble regimes
and large price volatility. The mapping of the detection of (financial)
instabilities to the measurement of the Kreiss constant improved conceptually
and operationally on the previous approaches attempting to anticipate critical
phase transitions [67, 13, 68, 69, 55], which do not incorporate the
ubiquitous non-normality of complex system dynamics.

## Conclusions

In most of the literature, financial and socio-economic bubbles have been
thought of as being associated with special regimes where self-reinforcing
interactions strengthen transiently towards a critical point and lead to some
form of collective exuberance. This has been formalized by models in physics,
ecology and mathematics assuming the presence of an underlying phase
transition, criticality, bifurcation or catastrophe. The focus put mainly on
critical points is particularly surprising insofar as the existence of non-
critical transients via non-normal dynamics has been well-established in other
fields of science such as in hydrodynamic turbulence.

Here, we have demonstrated, via agent-based model simulations and empirically,
that such transient phases of exuberance are also generic in real social
systems ubiquitously characterized by non-normal properties of asymmetric
hierarchical interactions and Ising-like imitative interactions. Our key
insight is that we find seemingly critical behavior in the sub-critical regime
as a result of transient collective behavior. This is reminiscent of the
behavior that would result in the agent-based model from the Ising-like
interactions between noise traders being at or close to criticality in a
normal network. In such normal networks of interacting noise traders, if the
Ising-like component is not crossing criticality, one does not observe
bubbles. It is the non-monotonic burst response to perturbation in the sub-
critical regime of the non-normal network of interactive noise traders that
create transient bubbles.

An important corollary is that financial bubbles should be expected as
intrinsic, rather than abnormalities appearing in very special conditions.
This may explain their ubiquity in financial markets, from past centuries to
the present. Our results have been derived based on first-order principles
that are not specific to our implementation but instead rely on the fact that,
on short time-scales, opinion-leaders can significantly lead opinions astray
from rational ones. With these general assumptions and the broad applicability
of models involving hierarchical and stochastic imitative interactions, our
framework is expected to explain ubiquitous crowd-forming patterns and
collective structures in general hierarchical social networks.

## Methods

### Dynamics of Noise Traders

In any ABM of financial markets, it is crucial to have at least two classes of
traders. When market participants become too homogeneous, liquidity rarifies
and no trades occur. We initiate the model with equal weight between the
fundamentalists and the noise traders, where their weights is determined by
their financial capital. As their wealth evolves with time, the relative
importance of each group varies with time. The specific description of the
investment strategy of fundamentalists is found in the Supplementary Note 1
and refs. [25, 26].

We label the set of $N$ noise traders (or nodes in the network) from
$1,\ldots,N$. We denote the state of each node $i=1,\ldots,N$ as
$s^{i}=\pm{1}$. The transition probability $\pi$ that trader $i$ flips its
state $s^{i}_{t}$ at time $t$ depends on the opinion of its in-neighbours,
according to

$\pi({s^{i}_{t+1}=-s^{i}_{t}})=\frac{p^{\pm}}{2}\left(1-\kappa\frac{1}{k^{in}_{i}}{s}^{i}_{t}\sum_{j}a_{ij}{s^{j}_{t}}\right)$
(4)

where $p^{\pm}$ controls the average holding time (net of social influence) of
each asset and the social coupling strength $\kappa$ determines the noise
traders’ susceptibility to social imitation. Here, $\mathbf{A}=\\{a_{ij}\\}$,
where $a_{ij}=1$ if there exists a directed edge (influence) from node $j$ to
node $i$, and $a_{ij}=0$ otherwise. The in-degree of node $i$ is the number of
directed edges pointing to node $i$, which is given by
$k_{i}^{in}=\sum_{j=1}^{N}a_{ij}$. In the Ising model, if node $i$ switches
its state from time step $t$ to time step $t+1$, i.e.
$s_{t+1}^{i}=-s_{t}^{i}$, the change of the value of node $i$’s state is
$-2s^{i}_{t}$. Given the probability of node $i$ to switch its state according
to (4), the average rate of change of the spin starting in the state
$s_{t}^{i}$ is given by $\Delta s^{i}_{t}=-2~{}s_{t}^{i}~{}\pi$. We introduce
the $N$-dimensional state column vector
$\vec{s}(t)=\left(s^{1}_{t},s^{2}_{t},...,s^{N}_{t}\right)$. Together with
(4), the average rate of state transition can then be written as
$\triangle{\vec{s}(t)}=\vec{s}(t+1)-\vec{s}(t)=p^{\pm}(\kappa\mathbf{\Lambda}\mathbf{A}-\mathbf{I})\vec{s}(t)$
where $\mathbf{\Lambda}$ is an $N\times{N}$ diagonal matrix with
$1\left/k^{in}_{i}\right.$ on the $i$-th diagonal entry and $\mathbf{I}$ is
the identity matrix. Introducing the matrix

$\mathbf{M}\equiv p^{\pm}(\kappa\mathbf{\Lambda}\mathbf{A}-\mathbf{I})~{},$
(5)

recovers the simple dynamical equation (1).

### A Primer on Non-Normality

Following a classic textbook [22], we briefly summarize the basic theory
behind non-normal matrices (see Supplementary Note 2 for details).

Let $\mathbf{M}$ be an $(N\times N)$-matrix. The set of all eigenvalues of
$\mathbf{M}$ is called the spectrum $\sigma(\mathbf{M})$. A matrix is called
normal if $\mathbf{M}^{T}\mathbf{M}=\mathbf{M}\mathbf{M}^{T}$, and the
spectral theorem asserts that each normal $\mathbf{M}$ has a set of $n$
pairwise orthonormal eigenvectors of $\mathbf{M}$. By contrast, if
$\mathbf{M}$ is non-normal,
$\mathbf{M}^{T}\mathbf{M}\neq\mathbf{M}\mathbf{M}^{T}$, there exist some
eigenvectors which are non-orthogonal to each other. Since symmetric matrices
are always normal, it is a necessary, but not a sufficient condition that
matrix (5) contains directed interactions to be considered non-normal.

If $\lambda$ is an eigenvalue of $\mathbf{M}$, the resolvent matrix
$\mathbf{M}-\lambda\mathbf{I}$ is not invertible since there exists an
eigenvector $\vec{v}$ with $(\mathbf{M}-\lambda\mathbf{I})\vec{v}=0$. An
alternative definition of the spectrum $\sigma(\mathbf{M})$ is thus the set of
points $\lambda\in\mathbb{C}$ where the resolvent matrix does not exit. But
the question “Does $\left(\mathbf{M}-\lambda\mathbf{I}\right)^{-1}$ exist?” is
binary and may change from “yes” to “no” by just a tiny
$\epsilon$-perturbation of $\lambda$. In the presence of noise, a better
question to ask is whether
$\left|\left|\left(\mathbf{M}-\lambda\mathbf{I}\right)^{-1}\right|\right|$ is
large with respect to some matrix norm $\left|\left|\cdot\right|\right|$. This
leads to the definition of the $\epsilon$-pseudospectrum, defined as the set
of points where
$\left|\left|\left(\mathbf{M}-\lambda\mathbf{I}\right)^{-1}\right|\right|$ is
large (larger then $\epsilon^{-1}$), or formally,
$\sigma_{\epsilon}(\mathbf{M})\equiv\left\\{\lambda\in\mathbb{C}~{}:~{}\left|\left|\left(\mathbf{M}-\lambda\mathbf{I}\right)^{-1}\right|\right|>\epsilon^{-1}\right\\}.$
The $\epsilon$-pseudospectrum is the open subset of the complex plane bounded
by the $\epsilon^{-1}$ level-curve of the norm of the resolvent. Intuitively,
one can then assume that the $\epsilon$-pseudospectrum is closely confined
around the eigenvalues of $\mathbf{M}$. For normal matrices, this assumption
is correct. However, for non-normal matrices it is not, and
$\left|\left|\left(\mathbf{M}-\lambda\mathbf{I}\right)^{-1}\right|\right|$ may
be large even when $\lambda$ is far away from the spectrum (Figure 1 (a) and
(b)).

Consider the proportional growth equation
$\left.\mathrm{d}\vec{s}\right/\mathrm{d}t=\mathbf{M}\vec{s}$ with explicit
solution $\vec{s}(t)=e^{t\mathbf{M}}\vec{s}(0)$. It is well-known that the
asymptotic behavior for $t\to\infty$ is governed by the largest real-part of
all eigenvalues of $\mathbf{M}$. For the short-term behavior, $t\downarrow 0$,
it can be shown that
$\frac{\mathrm{d}~{}}{\mathrm{d}t}\left|\left|e^{t\mathbf{M}}\right|\right|_{t=0}=\omega(\mathbf{M})\equiv\text{sup}~{}\sigma\left(\frac{1}{2}\left(\mathbf{M}+\mathbf{M}^{T}\right)\right)$
where $\omega(\mathbf{M})$ is called the numerical abscissa of $\mathbf{M}$
(Figure 1(c)). Our main interest are, however, intermediate values of $t$. To
describe such transient behavior, one has to consider the $\epsilon$-spectral
abscissa of a matrix $\mathbf{M}$ defined by
$\alpha_{\epsilon}(\mathbf{M})=\sup\text{Re}\left(\sigma_{\epsilon}(\mathbf{M})\right),$
i.e. the supremum of the real part of the $\epsilon$-pseudo-spectrum. An
important special case is the spectral abscissa
$\alpha(\mathbf{M})\equiv\alpha_{\epsilon=0}(\mathbf{M})$, defined as the
largest real-part of all eigenvalues of $\mathbf{M}$.

We now consider the case where $\alpha(\mathbf{M})<0$, i.e. where the long-
term behavior is asymptotically stable (Figure 1(c)), but
$\alpha_{\epsilon}(\mathbf{M})>0$ for some $\epsilon>0$. In that case, the
pseudospectra of $\mathbf{M}$ protrude significantly into the right-half plane
of $\mathbb{C}$, such that the real-parts of the pseudo-spectrum remain
positive (Figure 1 (e)). For any such non-normal matrix $\mathbf{M}$, the
Kreiss constant

$\mathcal{K}(\mathbf{M})\equiv\sup_{\epsilon>0}\frac{\alpha_{\epsilon}(\mathbf{M})}{\epsilon}$
(6)

is well-defined, and it can be shown that, for intermediate times (Figure
1(c)), there is transient growth according to (3).

In contrast to the ‘single-shot’ dynamics in which the system relaxes back to
the stable state after a short-lived shock, our model (1) is stochastic such
that the trajectory may get perturbed before it can equilibrate. It has been
shown that, in such stochastic models, non-normality will increase the
variance of the fluctuations, producing amplifications beyond what is expected
from its normal counterpart [70, 24, 42, 41].

Finally, we note that an alternative measure to characterize an operator’s
degree of non-normality is Henrici’s departure from normality

$d_{F}(\mathbf{A})=\left.\sqrt{\left|\left|\mathbf{A}\right|\right|_{F}^{2}-\sum_{\lambda\in\sigma(\mathbf{A})}\left|\lambda\right|^{2}}\right/\left|\left|\mathbf{A}\right|\right|_{F}^{2}.$
(7)

Henrici’s index is based on the observation that the Frobenius norm of a
normal matrix is given by
$\left|\left|\mathbf{A}\right|\right|_{F}^{2}=\text{tr}\left(\mathbf{A}^{T}\mathbf{A}\right)=\sum_{\lambda\in\sigma(\mathbf{A})}\left|\lambda\right|^{2}$.
It attains its minimum at zero once the matrix is normal and increases the
more the matrix deviates from normality. For example, the Henrici’s indices of
the adjacency matrices $\mathbf{A}$ of the networks depicted in Figure 2
(a,b,c) are equal to 1, 0.9007 and 0.8057, respectively.

As we show in the Supplementary Note 4 - at least within the realms of our
application - the relationship between $\mathcal{K}$ and $d_{F}$ is strictly
monotonous. Therefore all of our results presented below are also valid when
expressed in terms of $d_{F}(\mathbf{M})$.

### Parametrization of Non-Normal Matrices

The non-normal network with a total of $N$ nodes is initialized with $N_{0}$
top nodes. Each newly added node receives $m$ in-edges, i.e. channels of
communication through which it can be influenced. The source of each such edge
is selected with probability proportional to the existing nodes’ out-degree.
Once the $m$ source nodes of an added node are determined, each of the $m$
newly formed directed edges may be reciprocated with some probability $\rho$.
Based on empirical evidence from Reddit, we assume further that this
probability $\rho$ is modulated by the hierarchical level $\ell$ of each node,
$\rho=\rho(\ell)$. We thus model the rate of reciprocity as a sigmoid-function

$\rho(\ell)=\frac{\theta}{1+e^{-a(\ell(j)-b)}}~{}-~{}\underbrace{\frac{\theta}{1+e^{-a(1-b)}}}_{\equiv\gamma}$
(8)

where $\theta$ now serves as an upper asymptotic bound. In the remainder of
this article, we fix $\theta=0.2110,a=2.552$ and $b=3.668$ as determined
empirically below on meme-stock reddit data (Figure 2(d)). The offset
$\gamma\approx 0.0002$ has been added so that nodes in the highest level of
the hierarchy, $\ell=1$, are never reciprocated,
$\rho(\ell=1)=0~{}\forall\theta$. This offset is merely for convenience, so
that we can keep the number of top nodes (‘opinion leaders’) fixed. For large
values of $\ell$, we converge to the constant rate of reciprocity
$\theta-\gamma\approx\theta$ as in [48]. Parameter $\theta$ now has the
interpretation of the asymptotic level of reciprocity at high levels $\ell$. A
value of $\theta=1$ implies that most, albeit not all edges are reciprocated.
Details are found in Supplementary Note 3.

Our algorithm to generate non-normal adjacency matrix $\mathbf{A}$ has six
parameters: $N,N_{0},m,\theta,a$ and $b$. The price dynamics from (1) is not
directly governed by $\mathbf{A}$, but rather by the related matrix
$\mathbf{M}$. The two parameters $\kappa$ and $p^{\pm}$ allow us to control
the characteristics of $\mathbf{M}$ for given $\mathbf{A}$. In the remainder
of this article, we fix $p^{\pm}=0.05$ and we tune $\kappa$. We refer to the
Supplementary Note 4 for a detailed list of parameters, including parameters
for the fundamental agents and macro-economic quantities. These parameters
serve to fix the overall price scale, but do not qualitatively affect our
results. This leaves us with a three-parameter model, $N_{0},\theta$ and
$\kappa$. Importantly for what follows, the Kreiss constant $\mathcal{K}$ is
strictly decreasing in $\theta$ and increasing in $\kappa$, irrespective of
$N_{0}$, as long as $N_{0}\ll N$. Throughout this article, we constrain the
parameter such that $\alpha_{0}(\mathbf{M})<0$, i.e. the asymptotic system
dynamics is stable.

## References

  * [1] D. Sornette.  Critical phenomena in natural sciences (chaos, fractals, self-organization and disorder: Concepts and tools).  Springer Series in Synergetics, Heidelberg, 2004.
  * [2] R. Cont and J.-P. Bouchaud.  Herd behavior and aggregate fluctuations in financial markets.  Macroeconomic Dynamics, 4(2):170–196, 2000.
  * [3] M. Scheffer.  Critical transitions in nature and society.  Princeton University Press, 2009.
  * [4] J. C. Rocha, G. Peterson, Ö. Bodin, and S. Levin.  Cascading regime shifts within and across scales.  Science, 362(6421):1379–1383, 2018.
  * [5] T. M. Bury, R.I. Sujith, I. Pavithran, M. Scheffer, T. M. Lenton, M. Anand, and C. T. Bauch.  Deep learning for early warning signals of tipping points.  Proceedings of the National Academy of Sciences, 118(39), 2021.
  * [6] J. B. De Long, A. Shleifer, L. H. Summers, and R. J. Waldmann.  Positive feedback investment strategies and destabilizing rational speculation.  the Journal of Finance, 45(2):379–395, 1990.
  * [7] M. P. Scholl, A. Calinescu, and J. D. Farmer.  How market ecology explains market malfunction.  Proceedings of the National Academy of Sciences, 118(26), 2021.
  * [8] J.-P. Bouchaud and M. Potters.  Theory of financial risk and derivative pricing: from statistical physics to risk management.  Cambridge university press, 2003.
  * [9] D. Sornette.  Why stock markets crash: critical events in complex financial systems, volume 49.  Princeton University Press, 2nd printing, 2017.
  * [10] R. N. Mantegna and H. E. Stanley.  Introduction to econophysics: correlations and complexity in finance.  Cambridge University Press, 1999.
  * [11] S. Galam.  Sociophysics: A review of Galam models.  International Journal of Modern Physics C, 19(03):409–440, 2008\.
  * [12] S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley, and S. Havlin.  Catastrophic cascade of failures in interdependent networks.  Nature, 464(7291):1025–1028, 2010.
  * [13] M. Scheffer, S. R. Carpenter, T. M. Lenton, J. Bascompte, W. Brock, V. Dakos, J. Van de Koppel, I. A. Van de Leemput, S. A. Levin, E. H. Van Nes, et al.  Anticipating critical transitions.  Science, 338(6105):344–348, 2012.
  * [14] S. Battiston, J. D. Farmer, A. Flache, D. Garlaschelli, A. G. Haldane, H. Heesterbeek, C. Hommes, C. Jaeger, R. May, and M. Scheffer.  Complexity theory and financial regulation.  Science, 351(6275):818–819, 2016.
  * [15] D. Sornette.  Sweeping of an instability: an alternative to self-organized criticality to get power laws without parameter tuning.  J. Phys. I France, 4:209–221, 1994.
  * [16] B. Corominas-Murtra, J. Goñi, R. V. Solé, and C. Rodríguez-Caso.  On the origins of hierarchy in complex networks.  Proceedings of the National Academy of Sciences, 110(33):13316–13321, 2013.
  * [17] S.C. Lera and D. Sornette.  A theory of discrete hierarchies as optimal cost-adjusted productivity organisations.  PLoS ONE, 14(4):e0214911 (https://doi.org/10.1371/journal.pone.0214911), 2019.
  * [18] M. Asllani and T. Carletti.  Topological resilience in non-normal networked systems.  Physical Review E, 97(4):042302, 2018.
  * [19] J. D. O’Brien, K. A. Oliveira, J. P. Gleeson, and M. Asllani.  Hierarchical route to the emergence of leader nodes in real-world networks.  Physical Review Research, 3(2):023117, 2021.
  * [20] C. Duan, T. Nishikawa, D. Eroglu, and A. E. Motter.  Network structural origin of instabilities in large complex systems.  Science advances, 8(28):eabm8310, 2022.
  * [21] L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T. A. Driscoll.  Hydrodynamic stability without eigenvalues.  Science, 261(5121):578–584, 1993.
  * [22] M. Embree and L. N. Trefethen.  Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators.  Princeton University Press Princeton, 2005.
  * [23] B. K. Murphy and K. D. Miller.  Balanced amplification: a new mechanism of selective amplification of neural activity patterns.  Neuron, 61(4):635–648, 2009.
  * [24] T. Biancalani, F. Jafarpour, and N. Goldenfeld.  Giant amplification of noise in fluctuation-induced pattern formation.  Physical Review Letters, 118(1):018101, 2017.
  * [25] T. Kaizoji, M. Leiss, A. Saichev, and D. Sornette.  Super-exponential endogenous bubbles in an equilibrium model of fundamentalist and chartist traders.  Journal of Economic Behavior & Organization, 112:289–310, 2015\.
  * [26] R. Westphal and D. Sornette.  Market impact and performance of arbitrageurs of financial bubbles in an agent-based model.  Journal of Economic Behavior & Organization, 171:1–23, 2020.
  * [27] H. P. Minsky.  Monetary systems and accelerator models.  The American Economic Review, 47(6):860–883, 1957.
  * [28] C. Castellano, S. Fortunato, and V. Loreto.  Statistical physics of social dynamics.  Reviews of Modern Physics, 81(2):591, 2009.
  * [29] M. Gisler, D. Sornette, and R. Woodard.  Innovation as a social bubble: The example of the human genome project.  Research Policy, 40:1412–1425, 2011.
  * [30] D. Sornette.  Physics and financial economics (1776–2014): puzzles, Ising and agent-based models.  Reports on Progress in Physics, 77(6):062001, 2014.
  * [31] G. Dosi and A. Roventini.  More is different… and complex! the case for agent-based macroeconomics.  Journal of Evolutionary Economics, 29(1):1–37, 2019.
  * [32] T. Ott, P. Masset, T. S. Gouvea, and A. Kepecs.  Apparent sunk cost effect in rational agents.  Science Advances, 8(6):eabi7004, 2022.
  * [33] Didier Sornette and Peter Cauwels.  Financial bubbles: mechanisms and diagnostics.  Review of Behavioral Economics, 2(3):279–305, 2015.
  * [34] W. A. Brock and C. H. Hommes.  Heterogeneous beliefs and routes to chaos in a simple asset pricing model.  Journal of Economic dynamics and Control, 22(8-9):1235–1274, 1998\.
  * [35] W. A. Brock and S. N. Durlauf.  Discrete choice with social interactions.  The Review of Economic Studies, 68(2):235–260, 2001.
  * [36] C. Chiarella, R. Dieci, and X.-Z. He.  Heterogeneity, market mechanisms, and asset price dynamics.  In Handbook of financial markets: Dynamics and evolution, pages 277–344. Elsevier, 2009.
  * [37] M. G. Neubert and H. Caswell.  Alternatives to resilience for measuring the responses of ecological systems to perturbations.  Ecology, 78(3):653–665, 1997.
  * [38] N. Hatano and D. R. Nelson.  Localization transitions in non-hermitian quantum mechanics.  Physical Review Letters, 77(3):570, 1996.
  * [39] M. G. Neubert, H. Caswell, and J. D. Murray.  Transient dynamics and pattern formation: reactivity is necessary for turing instabilities.  Mathematical biosciences, 175(1):1–11, 2002.
  * [40] B. Ravoori, A. B. Cohen, J. Sun, A. E. Motter, T. E. Murphy, and R. Roy.  Robustness of optimal synchronization in real networks.  Physical review letters, 107(3):034102, 2011.
  * [41] S. Nicoletti, D. Fanelli, N. Zagli, M. Asllani, G. Battistelli, T. Carletti, L. Chisci, G. Innocenti, and R. Livi.  Resilience for stochastic systems interacting via a quasi-degenerate network.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(8):083123, 2019.
  * [42] S. Nicoletti, N. Zagli, D. Fanelli, R. Livi, T. Carletti, and G. Innocenti.  Non-normal amplification of stochastic quasicycles.  Physical Review E, 98(3):032214, 2018.
  * [43] R. Muolo, T. Carletti, J. P. Gleeson, and M. Asllani.  Synchronization dynamics in non-normal networks: the trade-off for optimality.  Entropy, 23(1):36, 2020.
  * [44] G. Hennequin, T. P. Vogels, and W. Gerstner.  Non-normal amplification in random balanced neuronal networks.  Physical Review E, 86(1):011909, 2012.
  * [45] E. Gudowska-Nowak, M. A. Nowak, D. R. Chialvo, J. K. Ochab, and W. Tarnowski.  From synaptic interactions to collective dynamics in random neuronal networks models: critical role of eigenvectors and transient behavior.  Neural Computation, 32(2):395–423, 2020.
  * [46] S. Tang and S. Allesina.  Reactivity and stability of large ecosystems.  Frontiers in Ecology and Evolution, 2:21, 2014.
  * [47] R. Muolo, M. Asllani, D. Fanelli, P. K. Maini, and T. Carletti.  Patterns of non-normality in networked systems.  Journal of theoretical biology, 480:81–91, 2019.
  * [48] M. Asllani, R. Lambiotte, and T. Carletti.  Structure and dynamical behavior of non-normal networks.  Science Advances, 4(12):eaau9403, 2018.
  * [49] G. Baggio, V. Rutten, G. Hennequin, and S. Zampieri.  Efficient communication over complex dynamical networks: The role of matrix non-normality.  Science Advances, 6(22):eaba2282, 2020.
  * [50] S. Johnson.  Digraphs are different: Why directionality matters in complex systems.  Journal of Physics: Complexity, 1(1):015003, 2020.
  * [51] M. Kawakatsu, P. S. Chodrow, N. Eikmeier, and D. B. Larremore.  Emergence of hierarchy in networked endorsement dynamics.  Proceedings of the National Academy of Sciences, 118(16), 2021.
  * [52] A. Fiske.  The four elementary forms of sociality: framework for a unified theory of social relations.  Psychological Review, 99(4):689, 1992.
  * [53] M. Favre and D. Sornette.  A generic model of dyadic social relationships.  PLoS ONE, 10(3):e0120882, 2015.
  * [54] A.-L. Barabási and R. Albert.  Emergence of scaling in random networks.  Science, 286(5439):509–512, 1999.
  * [55] S. C. Lera, A. Pentland, and D. Sornette.  Prediction and prevention of disproportionally dominant agents in complex networks.  Proceedings of the National Academy of Sciences, 117(44):27090–27095, 2020.
  * [56] S. Johnson and N. S. Jones.  Looplessness in networks is linked to trophic coherence.  Proceedings of the National Academy of Sciences, 114(22):5618–5623, 2017.
  * [57] C. Pilgrim, W. Guo, and S. Johnson.  Organisational social influence on directed hierarchical graphs, from tyranny to anarchy.  Scientific Reports, 4388(10), 2020.
  * [58] T. Lux.  Herd behaviour, bubbles and crashes.  The Economic Journal, 105(431):881–896, 1995.
  * [59] Zhi-Qiang Jiang, Wei-Xing Zhou, Didier Sornette, Ryan Woodard, Ken Bastiaensen, and Peter Cauwels.  Bubble diagnosis and prediction of the 2005–2007 and 2008–2009 chinese stock market bubbles.  Journal of Economic Behavior & Organization, 74(3):149–162, 2010\.
  * [60] A. Hüsler, D. Sornette, and C. H. Hommes.  Super-exponential bubbles in lab experiments: evidence for anchoring over-optimistic expectations on price.  Journal of Economic Behavior & Organization, 92:304–316, 2013\.
  * [61] T. M. Liggett.  Stochastic interacting systems: contact, voter and exclusion processes, volume 324.  springer science & Business Media, 1999.
  * [62] J. Miao and P. Wang.  Asset bubbles and credit constraints.  American Economic Review, 108(9):2590–2628, 2018.
  * [63] Š. Lyócsa, E. Baumöhl, and T. Vỳrost.  YOLO trading: Riding with the herd during the GameStop episode.  Finance Research Letters, page 102359, 2021.
  * [64] L. Lucchini, L. M. Aiello, L. Alessandretti, G. De Francisci Morales, M. Starnini, and Baronchelli A.  From Reddit to Wall Street: the role of committed minorities in financial collective action.  Royal Society Open Science, 9(211488), 2022.
  * [65] A. Betzer and J. P. Harries.  How online discussion board activity affects stock trading: the case of gamestop.  Financial Markets and Portfolio Management, pages 1–30, 2022.
  * [66] Z. G. Nicolaou, T. Nishikawa, S. B. Nicholson, J. R. Green, and A. E. Motter.  Non-normality and non-monotonic dynamics in complex reaction networks.  Physical Review Research, 2(4):043059, 2020.
  * [67] M. Scheffer, J. Bascompte, W. Brock, V. Brovkin, S. R. Carpenter, V. Dakos, H. Held, E. H. Van Nes, M. Rietkerk, and G. Sugihara.  Anticipating critical transitions.  Nature, 461:53–59, 2009.
  * [68] I. A. Van de Leemput, M. Wichers, A. O. J. Cramer, D. Borsboom, F. Tuerlinckx, P. Kuppens, Egbert H. van Nes, W. Viechtbauer, E. J. Giltay, S. H. Aggen, et al.  Critical slowing down as early warning for the onset and termination of depression.  Proceedings of the National Academy of Sciences, 111(1):87–92, 2014\.
  * [69] J. Ma, Y. Xu, Y. Li, R. Tian, and J. Kurths.  Predicting noise-induced critical transitions in bistable systems.  Chaos, 29(8):081102, 2019.
  * [70] B. F. Farrell and P. J. Ioannou.  Variance maintained by stochastic forcing of non-normal dynamical systems associated with linearly stable shear flows.  Physical Review Letters, 72:1188–1191, 1994.


# Models for the yielding behaviour of amorphous solids

Srikanth Sastry<EMAIL_ADDRESS>Jawaharlal Nehru Centre for Advanced
Scientific Research, Jakkur Campus, Bengaluru 560064, India.

(9-3-2020)

###### Abstract

Understanding the mechanical response and failure of solids is of obvious
importance in their use as structural materials. The nature of plastic
deformation leading to yielding of amorphous solids has been vigorously
pursued in recent years. Investigations employing both unidirectional and
cyclic deformation protocols reveal a strong dependence of yielding behaviour
on the degree of annealing. Below a threshold degree of annealing, the nature
of yielding changes qualitatively, to progressively more discontinuous
yielding. Theoretical investigations of yielding in amorphous solids have
almost exclusively focused on yielding under unidirectional deformation, but
cyclic deformation reveals several interesting features that remain largely
un-investigated. Focusing on athermal cyclic deformation, I investigate a
family of models based on an energy landscape description. These models
reproduce key interesting features observed in simulations, and provide an
interpretation for the intriguing presence of a threshold energy.

Amorphous solids of a wide variety are of scientific and technological
importance. Apart from hard glasses such as oxide glasses, metallic glasses,
etc., several soft materials, such as colloidal assemblies, gels, emulsions
and pastes also belong to the broad spectrum of amorphous solids Bonn  _et
al._ (2017); Nicolas  _et al._ (2018). The manner in which such solids respond
to external stresses is equally diverse, understanding of which is of interest
for many reasons – whereas the development of plastic deformation leading to
mechanical failure may be a phenomenon to predict and prevent in the case of
structural materials, controlling the elastoplastic flow properties is of
interest in the case of soft materials, such as gels and pastes Bonn  _et al._
(2017). A general description of the mechanical response of amorphous solids
needs to take note of the microscopic structural disorder and account for the
apparent diversity of responses.

The nature of elementary processes of plasticity, interactions between them,
the nature of the yielding transition and the dependence of the mechanical
response on preparation history, have been actively investigated over the
years Falk and Langer (2011); Dasgupta  _et al._ (2012); Lin and Wyart (2016);
Regev  _et al._ (2015); Jaiswal  _et al._ (2016); Jin  _et al._ (2018); Parisi
_et al._ (2017); Urbani and Zamponi (2017); Leishangthem  _et al._ (2017);
Ozawa  _et al._ (2018); Bhaumik  _et al._ (2019); Radhakrishnan and Fielding
(2016); Popović  _et al._ (2018); Barlow  _et al._ (2020). In particular,
computer simulations of athermal quasi-static shear (AQS) deformations Shi
_et al._ (2007); Leishangthem  _et al._ (2017); Ozawa  _et al._ (2018);
Bhaumik  _et al._ (2019); Yeh  _et al._ (2019); Vasisht  _et al._ (2020);
Priezjev (2018), employing both uniform and cyclic shear protocols, reveal
that the character of yielding depends strongly on the initial degree of
annealing of the solid. Under uniform shear, the evolution of stress with
strain is gradual for poorly annealing glasses, whereas well annealed glasses
display both a stress overshoot and a discontinuous jump down to the flow
stress, with the strength of discontinuity growing with annealing, below a
threshold. Such variation of yielding behaviour with annealing has been
suggested as a framework to understand the observed diversity among amorphous
solids Ozawa  _et al._ (2018).

In the case of cyclic shear, poorly annealed glasses, above a threshold degree
of annealing (or energy), anneal towards a threshold energy with increasing
amplitude of strain, before yielding Leishangthem  _et al._ (2017); Bhaumik
_et al._ (2019); Yeh  _et al._ (2019), accompanied by shear banding Parmar
_et al._ (2019) (as also for uniform shear). Initially well annealed samples
with energies below the threshold do not show any change in properties untill
they yield discontinuously. Above yielding, the properties do not depend on
the initial sample history. Interestingly, the number of cycles required to
reach the steady state increases strongly as the strain amplitude approaches,
and possibly diverges at, the yielding amplitudeFiocco  _et al._ (2013); Regev
_et al._ (2015); Leishangthem  _et al._ (2017). These features are
schematically summarised in Fig. 1.

Figure 1: (Top) Schematic representation of yielding behaviour under cyclic
shear. (a) (Left) Dependence of steady state energy on strain amplitude
$\gamma_{max}$ for different initial energies. (b) (Right) Time required to
reach the steady state vs. strain amplitude $\gamma_{max}$. (Bottom) A
schematic description of the model. The minimum ($E_{min}$) and maximum
($E_{max}$) energy, strain at which energy is minimum ($\gamma_{0}$) and the
stability range ($\Delta\gamma$) are indicated for one mesostate (black), from
which transitions are possible to two other mesostates shown (blue, orange)
but not to the third (red).

Based on the observations that plastic deformation involves spatially
localised arrangements (termed shear transformation zones), and that they lead
to the generation of long range stress fields, a variety of elasto-plastic
models Nicolas  _et al._ (2018) have been developed and studied, aiming to
capture macroscopic deformation response, through a combination of a coarse-
grained description of local plasticity, and a continuum description at larger
length scales. Attempts have been made to incorporate the role of annealing in
such models, for uniform shear Talamali  _et al._ (2012); Ozawa  _et al._
(2018); Popović  _et al._ (2018); Barlow  _et al._ (2020), but response to
cyclic shear remains largely unaddressed. In particular, features exhibited by
athermal cyclic shear, including mechanically induced aging or annealing, the
presence of a threshold energy across which the character of yielding changes,
the apparent divergence of time scales have not been demonstrated to arise in
such models. One may well demand that a satisfactory elasto-plastic model be
able to capture such phenomena. A desirable model should capture both athermal
and thermal effects, and have the potential to incorporate parametrically
properties that may be specific to a variety of amorphous solids. The ability
to capture memory effects seen in cyclically sheared glasses would be another
desirable feature Keim  _et al._ (2019); Mungan  _et al._ (2019). Such
expectations point to an energy based model Sollich  _et al._ (1997); Jagla
(2007); Merabia and Detcheverry (2016); Jagla (2017)) (such as the soft glassy
rheology or SGR model Sollich  _et al._ (1997)), rather than models defined in
terms of threshold stresses, as with typical elasto-plastic models.

The models studied here are defined in terms of the distribution of states or
energy minima that a mesoscale block of an amorphous solid can be in, and
rules for transitions between them. Inclusion of an extended array for such
blocks and interactions between them along previously developed lines Nicolas
_et al._ (2018) is an obvious next step but is not pursued here. Each state
(termed a mesostate and investigated in detail in Mungan  _et al._ (2019)) is
characterised by (i) an energy $E_{0}$ (taken to be negative always) at a
strain value $\gamma_{0}$ at which the energy is minimum, (ii) a stability
range in strain over which it is stable, and (iii) a form for the variation of
the energy with strain (see Fig. 1 for an illustration). These
characterisations can be made in principle in material specific ways.

For the present work, I choose the stability range to be
$\gamma_{\pm}=\gamma_{0}\pm\sqrt{-E_{0}}$, and the energy at a given strain
within the stability range as $E(\gamma,E_{0},\gamma_{0})=E_{0}+{\mu\over
2}(\gamma-\gamma_{0})^{2}$. These choices embody the expectation (supported by
numerical evidence Bhaumik  _et al._ (2019)) that the stability range as well
as the elastic energy increase before instability increase upon lowering the
minimum energy. I consider a Gaussian distribution of energies,
$P_{0}(E_{0})=\sqrt{2\over\pi\sigma^{2}}\exp(-{E_{0}^{2}\over 2\sigma^{2}})$,
for $-1<E_{0}<0$, $\sigma=0.1$ 111This choice, for concreteness, would
correspond to the density of states of a subvolume of $50$ particles, with the
number of states $\Omega\sim\exp(\alpha N)$, $\alpha\sim 1$ Sastry (2001).

For a given energy $E_{0}$, it is assumed that one has mesostates with several
possible values of the stress-free strain $\gamma_{0}$, reflecting the fact
that distinct configurations of the same energy may exist, related to each
other by a shift in the strain values at which they are stable. 222This in
principle introduces an additional $E_{0}$ dependent in the density of states,
but we consider an implementation where the total weight of all relevant
mesostates of a given energy is given by $P_{0}(E_{0})$.. For concreteness,
two specific cases are considered here: (1) The regular case, where mesostates
have $\gamma_{0}=n\times 2\sqrt{-E_{0}}$, $n$ being an integer, so that
mesostates at a given energy are present in a periodic, non-overlapping
fashion. (2) The uniform case, where the $\gamma_{0}$ values are uniformly
distributed.

Figure 2: (Left) The yielding diagram that shows the steady state energies
reached after repeated cyclic deformation, as a function of the strain
amplitude $\gamma_{max}$, for a range of initial energies. (Right) Along with
average energies over all samples (black), partial average energies are shown,
of samples that yield (red) and those that do not (blue), revealing the
discontinuous nature of the transition to the yielded state. The yielding
probability vs. strain is also shown (green).

When the stability limit $\gamma_{\pm}$ is reached, the state of the system
makes a transition from the current mesostate to another. In the athermal
case, clearly, such a transition is possible only to other states with lower
energy at the same strain value (see Fig. 1). Depending on the $E_{0}$ and
$\gamma_{0}$ values of the other states, a transition may occur to mesostates
with higher or lower $E_{0}$, a key feature on which the behaviour described
below depends (which requires $\mu>1$ for the regular case). One may consider
a transition rule that permits transitions to (i) any mesostate with a lower
energy at the transition strain, or as a physically motivated choice, (ii)
restrict the range of $E_{0}$ values to which transitions are permitted,
within a range $\delta E$. Both these cases are considered.

Figure 3: (Left Panel) Average energies as a function of the number of cycles.
(Right Panel) Relaxation time $\tau$ to reach the steady state below yielding.
The dashed line is a power law fit, and the solid line describes
$\tau=A~{}\gamma_{max}\exp({\gamma_{max}^{4}/2\sigma^{2}})$. The inset shows a
log-log plot of $\tau$ vs. $\gamma_{max}\exp({\gamma_{max}^{4}/2\sigma^{2}})$.

The case of regular distribution of $\gamma_{0}$ is considered first, with a
choice of $\mu=1.1$ and investigated numerically as follows. Starting with an
initial energy $E_{0}$, and $\gamma_{0}=0$, the strain $\gamma$ is varied
cyclically with amplitude $\gamma_{max}$, and the state of the system is
propagated untill a stability limit $\gamma_{\pm}$ is reached. When this
happens, the set of all states to which a transition can occur, with the
corresponding weights (given by $P_{0}(E_{0})$; in practice, the $E_{0}$ range
is divided into $10^{4}$ bins) is evaluated, and a new state, with its $E_{0}$
and $\gamma_{0}$ values, is randomly chosen among permitted states. The
procedure is repeated for $10^{3}$ cycles of shear. If the state of the system
does not change at the end of successive cycles, the procedure is terminated.
Otherwise, the zero strain energies for the last $200$ cycles are averaged.
This procedure is repeated for $10^{3}$ independent sample runs. Fig. 2 (a)
shows the resulting final energies as a function of both a set of initial
energies and strain amplitudes. Strikingly, all the essential features of
observations from simulations are reproduced. For high initial energy values,
energies shift to lower energies, towards a threshold value common to a range
of initial energies, before a discontinuous transition occurs, after which the
energies follow a common branch regardless of initial energy values. For
initial energies lower than the threshold, there is no change vs.
$\gamma_{max}$ till a yield value is reached, at which a discontinuous jump in
energy is seen. Naturally, the discontinuity in energy grows with the
reduction of the initial energy. From Fig. 2 (a), the discontinuous nature of
the transition is not fully apparent. To reveal it more clearly, for any
choice of $\gamma_{max}$, the final energies for the cases where a transition
occurs (and those for which it does not) are computed, as also the probability
of occurrence of the transition vs. $\gamma_{max}$. These are shown for
$E_{0}^{init}=-0.15$ in Fig. 2 (b), which illustrates that (i) the yielding
probability rises sharply from $0$ to $1$ in a very narrow range of
$\gamma_{max}$, and (ii) the energy averaged over cases where a transition
does not occur continues to decrease with $\gamma_{max}$ even beyond the yield
strain range, whereas when a transition occurs, a nearly constant final energy
is reached. The latter observation, interestingly, echoes simulation results
Parmar  _et al._ (2019) where continued annealing is observed beyond the yield
strain amplitude, away from the shear band in which strain gets localised. To
interpret the results here suitably, one must consider that the present system
represents the behaviour of one mesoscale block in a more extended system.
Better still, such an extended system should be investigated, which will be
taken up in future work. A comment on the average energies beyond yielding is
also in order in this context. From the results here, the post-yield energy
does not continue to increase, whereas in simulations, it does. This can be
simply interpreted as a result of a larger fraction of the simulated system
undergoing yielding rather than a change in the typical energy within that
subvolume. Results in Parmar  _et al._ (2019) clearly support such an
interpretation.

Next, the energies of individual samples vs. cycle number averaged for
$\gamma_{max}$ below the yield value, are shown in Fig. 3 for an initial
energy of $-0.15$. The time scales (measured in terms of cycles) to reach the
final state increase with $\gamma_{max}$ consistently with simulation results
Fiocco  _et al._ (2013); Leishangthem  _et al._ (2017), and can be fitted to a
power law as shown, but the fit parameters are not meaningful. A way of
understanding the increasing time scales is presented below, which leads to
the expectation that $\tau\sim\gamma_{max}\exp({\gamma_{max}^{4}\over
2\sigma^{2}})$ which provides a very reasonable description. No long
relaxation times are observed above the transition, which I interpret to
indicating that the long relaxation times observed in simulations arise from
parts out the system outside the yielded region rather than the yield region
itself, which can be tested in simulations.

Figure 4: Yielding diagram for different cases: (a) The constrained case with
$\delta E=0.05$, $\sigma=0.1$, $\mu=1.1$ . (b) The uniform case with
$\sigma=0.1$, $\mu=1.1$ (c) The regular case with $\sigma=0.15$, $\mu=1.1$ (d)
The SGR version, with $\sigma=0.1$, $\mu=2$.

Before discussing the results further, I consider variations of the model
described above. The first is the regular model discussed above, but with a
constraint on the range of energy change when a mesostate transition occurs
(with the choice $\delta E=0.05$, $\sigma=0.1$, $\mu=1.1$). The second is the
uniform model ($\sigma=0.1$, $\mu=1.1$) with $\gamma_{0}$ values uniformly
distributed, rather than discretely, for any $E_{0}$. I also consider the
regular case with $\sigma=0.15$, $\mu=1.1$. Finally, with $\mu=2$, which
permits transitions to all mesostates without regard to their energy (since
the threshold energy for all states is $0$), and with the additional condition
that the $\gamma_{0}$ value for the new state is the same as the current value
of strain, one realises a specific instance that is the same as the SGR model
Sollich  _et al._ (1997). The relevant yielding diagrams are reported Fig. 4.
Investigation of these variations reveals that qualitatively, each of them
reproduces the behaviour discussed above. Equally importantly, these
variations do change the threshold energy as well as the energy beyond
yielding, which may be relevant in understanding material specific behaviour.
The constrained case demonstrates that transitions to widely different states
is not a necessary feature to produce the observed behaviour. The regular case
with $\sigma=0.15$ illustrates that the distribution $P_{0}(E_{0})$ has a
strong bearing on the threshold energy. Importantly, it also suggests that
models with a trivial (e. g., flat) distributions in the stability range are
likely not to capture the yielding transition under athermal cyclic shear,
which should be tested. Indeed, for $\sigma=0.25$, the simulations converge to
energies below $-\gamma_{max}^{2}$ all the way to the lowest energies
considered. The results indicate that in the present formulation, limit cycles
associated with memory effects Adhikari and Sastry (2018); Mungan  _et al._
(2019) do not arise. Constructing suitable transition graphs Mungan  _et al._
(2019); Mungan and Witten (2019); Fiocco  _et al._ (2015) that can describe
limit cycles associated with memory effects need to be explored.

Figure 5: The histogram of energies reached after one cycle (y-axis) vs. the
initial energy (x-axis) for a range of strain amplitudes $\gamma_{max}$
indicated in the panels, for the uniform case with $\mu=1.1$. The yield
$\gamma_{max}$ in this case is $\approx 0.563$.

The most intriguing result here and in earlier simulations is the presence of
the threshold energy, which is not an obvious part of the model description.
In contrast, the mean energy beyond yielding has the simple explanation. It is
the energy above which the mean energy to which a mesostate transitions is
lower, whereas below, it is higher. To investigate the meaning of the
threshold energy, I consider, given a strain amplitude $\gamma_{max}$, and
initial energy, the distribution of energies to which the system transitions,
at the end of a full cycle, for the uniform case with $\mu=1.1$. Such two
dimensional distributions $P(E_{0},E_{1})$ are shown for a series of
$\gamma_{max}$ values in Fig. 5. For any $\gamma_{max}$, there is trivially a
range of energies for which no transition takes place, but for energies above,
the transition probabilities are roughly independent of the initial energy
(though this feature is not common to all the cases studied, behaviour in
other cases is qualitatively the same). It can be seen that the probability to
transition to an energy below the threshold for the given strain
($=-\gamma_{max}^{2}$, indicated by a horizontal line) decreases dramatically
with an increase in $\gamma_{max}$, and, between $\gamma_{max}=0.5$ and $0.7$,
it drops strongly (the strain at which the yielding probability is $1/2$ in
this case is $\approx 0.563$). Thus, the probability to transition to the
yielded state becomes overwhelming, but is never equal to $1$.
Correspondingly, the probability over a large number of cycles to transition
to the yielded state becomes a sharper function, but moves to higher strain
values.

Keeping in mind the case when $\mu=2$ wherein an unstable mesostate can make a
transition to any other mesostate without restriction, a simple, if
approximate estimate of the transition probability to yielding can be made as
follows. Since the stability range for a state of energy $-E_{0}$ is
$2\sqrt{-E_{0}}$, over a cycle of strain (or any unit of strain), the
transition rate out of that state is inversely proportional to
$\sqrt{-E_{0}}$, and the probability to transition to a new state of energy
$E^{{}^{\prime}}$ is $P_{0}(E^{{}^{\prime}})$. If one assumes that as a result
of these transitions, a stationary distribution of reached within a single
cycle, the stationary distribution is given by
$P_{st}(E_{0})\propto\sqrt{-E_{0}}P(E_{0})$. This is inaccurate for low
energies, and also does not incorporate the fact that once an energy below
$-\gamma_{max}^{2}$ is reached, no further transitions occur. Nevertheless, it
is a good approximate description of the distribution of energies reached
after one cycle. One can evaluate the probability, after one cycle, of getting
trapped in a minimum that is stable at $\gamma_{max}$,
$p_{trap}(\gamma_{max})=\int_{-\infty}^{-\gamma_{max}^{2}}P_{st}(E_{0})dE_{0}=\Gamma({3\over
4},{\gamma_{max}^{4}\over 2\sigma^{2}})/\Gamma({3\over
4})\sim\exp(-{\gamma_{max}^{4}\over 2\sigma^{2}})/\gamma_{max}$. The
relaxation time should go inversely as $p_{trap}$,
$\tau\sim\gamma_{max}\exp({\gamma_{max}^{4}\over 2\sigma^{2}})$, which was
seen to be good description of the simulated results (Fig. 3 (b)). Writing the
probability after $N_{cyc}$ cycles to yield as
$p_{trans}=(1-p_{trap})^{N_{cyc}}\sim\exp(-p_{trap}N_{cyc})$ we see that the
$p_{trans}$ will be an increasingly sharp function of $\gamma_{max}$ and the
yield strain amplitude, $\gamma_{max}^{Y}\sim(\log(N_{cyc}))^{1/4}$, depends
very weakly on $N_{cyc}$. Whether interactions between mesoscale regions in a
macroscopic sample could lead to a sharp transition needs to be investigated.
In Bhaumik  _et al._ (2019), it was observed that the threshold energy closely
corresponds to the energy at which the average index of saddle points goes to
zero, marking a change in the character of the energy landscape. The
implications of such a qualitative change, and how it may be incorporated in
the present description, merits further investigation.

In summary, the yielding behaviour under athermal cyclic deformation has been
investigated for a family of models envisaged to represent a mesoscopic
system, which remarkably captures several key observations for model glasses
that have been investigated through computer simulations using athermal
quasistatic shear. The mechanism for the transition is a competition of the
stability against cyclic shear of low energy mesostates, and the entropic
drive associated with the large number of higher energy minima present,
reminiscent of arguments presented in Jaiswal  _et al._ (2016); Parisi  _et
al._ (2017). The transition is always discontinuous, as also noted for uniform
shear in a recent theoretical analysis Barlow  _et al._ (2020).These suggest
that such a description, embedded in an elasto-plastic scheme is suitable for
further investigation as an approach to modelling the larger length scale
behaviour of sheared amorphous solids, such as the emergence of permanent or
multiple shear bands Vasisht  _et al._ (2018). As the approach is based on a
local energy landscape description, it has advantages in capturing thermal
behaviour as well as athermal behaviour. It also suggests detailed
investigation of the energetic aspects of local plastic events in simulations
Liu  _et al._ (2020) in order to understand better how a diversity of
amorphous solids may be satisfactorily modelled.

###### Acknowledgements.

I thank Vishwas Vasisht, Muhittin Mungan, Pallabi Das, Monoj Adhikari and
Surajit Sengupta for comments on the manuscript, Yagyik Goswami with help with
manuscript preparation, and support through the J. C. Bose Fellowship, SERB,
DST, India.

## References

  * Bonn  _et al._ (2017) D. Bonn, M. M. Denn, L. Berthier, T. Divoux, and S. Manneville, Reviews of Modern Physics 89, 035005 (2017).
  * Nicolas  _et al._ (2018) A. Nicolas, E. E. Ferrero, K. Martens, and J.-L. Barrat, Reviews of Modern Physics 90, 045006 (2018).
  * Falk and Langer (2011) M. L. Falk and J. S. Langer, Annu. Rev. Condens. Matter Phys. 2, 353 (2011).
  * Dasgupta  _et al._ (2012) R. Dasgupta, H. G. E. Hentschel, and I. Procaccia, Physical Review Letters 109, 255502 (2012).
  * Lin and Wyart (2016) J. Lin and M. Wyart, Physical review X 6, 011005 (2016).
  * Regev  _et al._ (2015) I. Regev, J. Weber, C. Reichhardt, K. A. Dahmen, and T. Lookman, Nature Communications 6, 8805 (2015).
  * Jaiswal  _et al._ (2016) P. K. Jaiswal, I. Procaccia, C. Rainone, and M. Singh, Physical Review Letters 116, 85501 (2016).
  * Jin  _et al._ (2018) Y. Jin, P. Urbani, F. Zamponi, and H. Yoshino, Science Advances 4 (2018).
  * Parisi  _et al._ (2017) G. Parisi, I. Procaccia, C. Rainone, and M. Singh, Proceedings of the National Academy of Sciences 114, 5577 (2017).
  * Urbani and Zamponi (2017) P. Urbani and F. Zamponi, Physical Review Letters 118, 38001 (2017).
  * Leishangthem  _et al._ (2017) P. Leishangthem, A. D. S. Parmar, and S. Sastry, Nature Communications 8, 14653 (2017).
  * Ozawa  _et al._ (2018) M. Ozawa, L. Berthier, G. Biroli, A. Rosso, and G. Tarjus, Proceedings of the National Academy of Sciences 115, 6656 (2018).
  * Bhaumik  _et al._ (2019) H. Bhaumik, G. Foffi, and S. Sastry, arXiv preprint arXiv:1911.12957 (2019).
  * Radhakrishnan and Fielding (2016) R. Radhakrishnan and S. M. Fielding, Physical Review Letters 117 (2016).
  * Popović  _et al._ (2018) M. Popović, T. W. de Geus, and M. Wyart, Physical Review E 98, 040901 (2018).
  * Barlow  _et al._ (2020) H. J. Barlow, J. O. Cochran, and S. M. Fielding, Phys. Rev. Lett. 125, 168003 (2020).
  * Shi  _et al._ (2007) Y. Shi, M. B. Katz, H. Li, and M. L. Falk, Physical review letters 98, 185505 (2007).
  * Yeh  _et al._ (2019) W.-T. Yeh, M. Ozawa, K. Miyazaki, and L. Berthier, arXiv preprint arXiv:1911.12951 (2019).
  * Vasisht  _et al._ (2020) V. V. Vasisht, G. Roberts, and E. Del Gado, Phys. Rev. E 102, 010604 (2020).
  * Priezjev (2018) N. V. Priezjev, Computational Materials Science 150, 162 (2018).
  * Parmar  _et al._ (2019) A. D. Parmar, S. Kumar, and S. Sastry, Physical Review X 9, 021018 (2019).
  * Fiocco  _et al._ (2013) D. Fiocco, G. Foffi, and S. Sastry, Physical Review E 88, 020301 (2013).
  * Talamali  _et al._ (2012) M. Talamali, V. Petäjä, D. Vandembroucq, and S. Roux, Comptes Rendus Mecanique 340, 275 (2012).
  * Keim  _et al._ (2019) N. C. Keim, J. D. Paulsen, Z. Zeravcic, S. Sastry, and S. R. Nagel, Rev. Mod. Phys. 91, 035002 (2019).
  * Mungan  _et al._ (2019) M. Mungan, S. Sastry, K. Dahmen, and I. Regev, Phys. Rev. Lett. 123, 178002 (2019).
  * Sollich  _et al._ (1997) P. Sollich, F. Lequeux, P. Hébraud, and M. E. Cates, Physical review letters 78, 2020 (1997).
  * Jagla (2007) E. Jagla, Physical Review E 76, 046119 (2007).
  * Merabia and Detcheverry (2016) S. Merabia and F. Detcheverry, EPL (Europhysics Letters) 116, 46003 (2016).
  * Jagla (2017) E. Jagla, Physical Review E 96, 023006 (2017).
  * Note (1) This choice, for concreteness, would correspond to the density of states of a subvolume of $50$ particles, with the number of states $\Omega\sim\mathop{exp}\nolimits(\alpha N)$, $\alpha\sim 1$ Sastry (2001).
  * Note (2) This in principle introduces an additional $E_{0}$ dependent in the density of states, but we consider an implementation where the total weight of all relevant mesostates of a given energy is given by $P_{0}(E_{0})$.
  * Adhikari and Sastry (2018) M. Adhikari and S. Sastry, The European Physical Journal E 41, 105 (2018).
  * Mungan and Witten (2019) M. Mungan and T. A. Witten, Phys. Rev. E 99, 052132 (2019).
  * Fiocco  _et al._ (2015) D. Fiocco, G. Foffi, and S. Sastry, Journal of Physics: Condensed Matter 27, 194130 (2015).
  * Vasisht  _et al._ (2018) V. V. Vasisht, M. L. Goff, K. Martens, and J.-L. Barrat, arXiv preprint arXiv:1812.03948 (2018).
  * Liu  _et al._ (2020) C. Liu, S. Dutta, P. Chaudhuri, and K. Martens, arXiv preprint arXiv:2007.07162 (2020).
  * Sastry (2001) S. Sastry, Nature 409, 164 (2001).