# Ordering dynamics and aging in the Symmetrical Threshold model

David Abella<EMAIL_ADDRESS>Instituto de Física Interdisciplinar y
Sistemas Complejos IFISC (CSIC - UIB), Campus UIB, 07122 Palma de Mallorca,
Spain  Juan Carlos González-Avella  Instituto de Física Interdisciplinar y
Sistemas Complejos IFISC (CSIC - UIB), Campus UIB, 07122 Palma de Mallorca,
Spain  Advanced Programming Solutions SL, 07121 Palma de Mallorca, Spain  Maxi
San Miguel  Instituto de Física Interdisciplinar y Sistemas Complejos IFISC
(CSIC - UIB), Campus UIB, 07122 Palma de Mallorca, Spain  José J. Ramasco
Instituto de Física Interdisciplinar y Sistemas Complejos IFISC (CSIC - UIB),
Campus UIB, 07122 Palma de Mallorca, Spain

###### Abstract

The so-called Granovetter-Watts model was introduced to capture a situation in
which the adoption of new ideas or technologies requires a certain redundancy
in the social environment of each agent to take effect. This model has become
a paradigm for complex contagion. Here we investigate a symmetric version of
the model: agents may be in two states that can spread equally through the
system via complex contagion. We find three possible phases: a mixed one
(dynamically active disordered state), an ordered one, and a heterogeneous
frozen phase. These phases exist for several configurations of the contact
network. Then we consider the effect of introducing aging as a non-Markovian
mechanism in the model, where agents become increasingly resistant to change
their state the longer they remain in it. We show that when aging is present,
the mixed phase is replaced, for sparse networks, by a new phase with
different dynamical properties. This new phase is characterized by an initial
disordering stage followed by a slow ordering process towards a fully ordered
absorbing state. In the ordered phase, aging modifies the dynamical
properties. For random contact networks, we develop a theoretical description
based on an Approximate Master Equation that describes with good accuracy the
results of numerical simulations for the model with and without aging.

††preprint: APS/123-QED

##  I  Introduction

A variety of collective phenomena can be well understood through stochastic
binary-state models for interacting agents. In these models, each agent is
assumed to be in one of two possible states, such as susceptible/infected,
adopters/non-adopters, etc., depending on the context of the model. The
interaction among agents is determined by the underlying contact network and
the dynamical rules of the model. There are various examples of binary-state
models, including processes of opinion formation [1, 2, 3, 4, 5] and disease
or social contagion [6, 7], among others. The consensus problem consists in
determining under which circumstances the agents end up sharing the same state
or when the coexistence of both states prevails. This is characterized by a
phase diagram that provides the boundaries separating domains of different
behaviors in the control parameter space. Macroscopic descriptions of these
models in terms of mean-field, pair, and higher-order approximations are well
established [8].

An important category of binary-state models are threshold models [9], which
were originally introduced by M. Granovetter [6] to address problems of social
contagion such as rumor propagation, innovation adoption, riot participation,
etc. Multiple exposures, or group interaction, are necessary in these models
to update the current state, a characteristic of complex contagion models [10,
11]. The threshold model presents a discontinuous phase transition from a
“global cascade” phase to a “no cascade” phase, which was analyzed in detail
in Ref. [9]. This model has been extensively studied on various network
topologies, such as regular lattices, small-world [10], random [12], clustered
[13, 14], modular [15], hypergraphs [16], homophilic [17] and coevolving [18]
networks.

A main difference between the threshold model and other binary-state models,
such as the Voter [1], majority vote (MV) [19, 20, 21], and nonlinear Voter
model [22, 23, 24, 25, 26, 27], is the lack of symmetry between the two
states. In the threshold model, changing state is only possible in one
direction, representing the adoption forever of a new state that initially
starts in a small minority of agents. A symmetric version of the threshold
model, with possible changes of states in both directions, was introduced in
Refs. [28, 29] to investigate the impact of noise and anticonformity. However,
a complete characterization of the Symmetrical Threshold model and its
ordering dynamics have not been addressed so far.

Aging is an important non-Markovian effect in binary-state models that has
significant implications. It describes how the persistence time of an agent in
a particular state influences the transition rate to a different state [30,
31, 32, 33, 34]. As such, the longer an agent remains in the current state,
the smaller the probability of changing. Aging has been shown to cause
coarsening dynamics towards a consensus state in the Voter model [31, 35], to
induce bona fide continuous phase transitions in the noisy Voter model [36,
37], modify the phase diagram and non-equilibrium dynamics of the Schelling
segregation model [38], and to modify non-trivially the cascade dynamics of
the threshold model [39]. The introduction of aging is motivated by strong
empirical evidence that human interactions do not occur at a constant rate and
cannot be described using a Markovian assumption. Empirical studies have
reported heavy-tail inter-event time distributions that reflect heterogeneous
temporal activity patterns in social interactions [40, 41, 42, 43, 44, 45].

In this work, we present a comprehensive analysis of the Symmetrical Threshold
model, including its full phase diagram, and we investigate the effects of
aging in the model. The model is examined in various network topologies, such
as the complete graph, Erdős-Rényi (ER) [46], random regular (RR) [47], and a
two-dimensional Moore lattice. The possible phases of the system are defined
by the final stationary state as well as by the ordering/disordering dynamics
characterized by the time-dependent magnetization and interface density, the
persistence, and the mean internal time. For both the original model and the
aging variant, the results of Monte Carlo numerical simulations are compared
with results from the theoretical framework provided by an Approximate Master
Equation (AME)[48, 39], which is general for any random network. We also
derive a mean-field analysis to describe the outcomes in a complete graph.

The article is organized as follows: In Section II, we describe the
Symmetrical Threshold model and provide the numerical and theoretical analysis
of the phase diagram. Each subsection reports the results for the different
networks chosen. Section III presents the Symmetrical Threshold model with
aging, the corresponding numerical and theoretical analysis, and the
comparison with the model without aging. The results for the Moore lattice are
shown in Section IV. Finally, we conclude with a summary and conclusions in
Section V.

##  II  Symmetrical Threshold model

Figure 1:  (a) Phase diagram of the Symmetrical Threshold model in a Complete
graph of $N=2500$ nodes. Dotted and dashed lines correspond to
$T=(1-|m_{0}|)/2$ and $T=(1+|m_{0}|)/2$, respectively. Average performed over
5000 realizations. (b) Potential representation from Eq. (II.1) for a set of
values of the threshold $T$, shown in different colors.

The system consists of a set of $N$ agents located at the nodes of a network.
The variable describing the state of each agent $i$ takes one of the two
possible values: $s_{i}=\pm 1$. Every agent has assigned a fixed threshold
$0\leq T\leq 1$, which determines the fraction of different neighbors required
to change state. Even though this value might be agent dependent, we will
consider here only the case with a homogeneous $T$ value for all the agents of
the system. In each update attempt, an agent $i$ (called active agent) is
randomly selected, and if the fraction of neighbors with a different state is
larger than the threshold $T$, the active agent changes state $s_{i}\to-
s_{i}$. Simulation time is measured in Monte Carlo (MC) steps, i.e., $N$
update attempts. Numerical simulations run until the system reaches a frozen
configuration (absorbing state) or until the average magnetization,
$m=(1/N)\sum_{i}s_{i}$, fluctuates around a constant value.

###  II.1 Mean-field

We first consider the mean-field case of the complete graph (all-to-all
connections). We take an initial random configuration with magnetization
$m_{0}$ and perform numerical simulations for various values of $T$ to
construct the phase diagram (shown in Fig. 1a). We find three different phases
based on the final state:

  * •

Phase I or Mixed: The system reaches an active disordered state (final
magnetization $m_{f}=0$) where the agents change their state continuously;

  * •

Phase II or Ordered: The system reaches the ordered absorbing states
($m_{f}=\pm 1$) according to the initial magnetization $m_{0}$;

  * •

Phase III or Frozen: The system freezes at the initial random state
$m_{f}=m_{0}$.

Figure 2:  Phase diagram of the Symmetrical Threshold model in an ER (a) and a
RR (b) graph, both of $N=4\cdot 10^{4}$ nodes and mean degree $\langle
k\rangle=8$. The color map indicates the value of the average final
magnetization $m_{f}$. The red dashed line is the HMF prediction of the mixed-
ordered critical line. The black solid lines correspond to the AME prediction
of the borders of Phase II. (c) Average final magnetization $m_{f}$ as a
function of the initial magnetization $m_{0}$ for different $T$ values
(indicated with different colors and markers) in the RR graph. The average is
performed over 5000 realizations. The dotted and solid lines are the HMF (for
$T=1/8-4/8$) and AME predictions (for all $T$), respectively.

For a given initial magnetization $m_{0}\neq 0$ and increasing $T$, the system
undergoes a mixed-ordered transition at a critical threshold
$T_{c}=(1-|m_{0}|)/2$, and an ordered-frozen transition at a critical
threshold $T_{c}^{*}=(1+|m_{0}|)/2>T_{c}$ (indicated by dotted and dashed
black lines in Fig. 1a). In this mean-field scheme, if the fraction of nodes
in state $+1$ is denoted by $x$, the condition for a node in state $-1$ to
change its state is given by $\theta(x-T)$, where $\theta$ is the Heaviside
step function. Thus, in the thermodynamic limit ($N\to\infty$), the variable
$x$ evolves over time according to the following mean-field equation:

$\frac{dx}{dt}=(1-x)\;\theta(x-T)-x\;\theta(1-x-T)=-\frac{\partial
V(x)}{\partial x}.$ (1)

Here, $V(x)$ is the potential function. The stationary value of $x$, $x_{\rm
st}$, is the solution of the implicit equation resulting from setting the time
derivative equal to 0. No analytical expression of $x_{\rm st}$ has been found
in terms of $T$, but the solutions can be understood in terms of the potential
$V(x)$:

$\displaystyle V(x)=$
$\displaystyle-\int(1-x)\;\theta(x-T)-x\;\theta(1-x-T)\;dx$ $\displaystyle=$
$\displaystyle\frac{x^{2}}{2}+\frac{1}{2}\left(T^{2}-2T-x^{2}+1\right)\;\theta(T+x-1)$
$\displaystyle-\frac{1}{2}\left(T^{2}-2T-x(x-2)\right)\;\theta(x-T).$ (2)

The minimum and maximum values of $V(x)$ correspond to stable and unstable
solutions, respectively. Figure 1b shows the potential’s dependence on the
magnetization, obtained after a variable change $m=2x-1$ in Eq. (II.1). For
$T<0.5$, $m=0$ is a stable solution, but increasing the threshold reduces the
range of values of the initial magnetization from which this solution is
reached, enclosing Phase I between the unstable solutions $m=1-2\,T$ and
$2\,T-1$. In fact, if $m_{0}>1-2\,T$, the system reaches the absorbing
solution $m=+1$, while if $m_{0}<-1+2\,T$, it reaches $m=-1$ (Phase II). For
$T=0.5$, there is just one unstable solution at $m=0$, and all the initial
magnetization values reach the absorbing states $m=\pm 1$. For $T>0.5$, the
potential is equal to a constant value for a range of $m_{0}$, which means
that an initial condition will remain in this state forever (Phase III). The
range of values of the initial condition from which this phase is reached
grows linearly with $T$ until $T=1$, where all initial conditions fulfill
$\frac{dm}{dt}=0$.

Note that the mean-field Symmetrical Threshold model for $T=1$ shows the same
potential profile as the mean-field Voter model [3, 1, 49]. The important
difference is that for the Voter model, any initial magnetization is
marginally stable, while in our model any initial magnetization is an
absorbing state in Phase III. In the Voter model finite size fluctuations will
take the system to the absorbing states $m=\pm 1$.

###  II.2 Random networks

We analyze the phase diagram of the Symmetrical Threshold model in two random
networks: Erdős-Rényi (ER) [46] and random regular (RR) [47] graphs with mean
degree $\langle k\rangle=8$. Figures 2a and 2b show the phase diagram for both
networks, where it is shown that the existence of the three phases previously
described is robust to changes in network structure. The main difference from
the all-to-all scenario is that Phase III does not freeze exactly at the same
initial magnetization. Instead, the system reaches an absorbing state with a
higher magnetization $m_{f}>m_{0}$. In this phase, the value of $m_{f}$
depends on the threshold such that increasing $T$, increases the disorder in
the system, until $T=1$, where $m_{f}=m_{0}$ (see Fig. 2c). On the other hand,
phases I and II reach the same stationary state as in the mean-field case.
Furthermore, the critical thresholds $T_{c}$ and $T_{c}^{*}$ show a different
dependence on $m_{0}$ depending on the network structure.

To explain the transitions exhibited by the model, we use a theoretical
framework for binary-state dynamics in complex networks [48]: the Approximate
Master Equation (AME), which considers agents in both states $\pm 1$ with
degree $k$, $m$ neighbors in state $-1$ that have been $j$ time steps in the
current state (called “internal time” or “age”) as different sets in a
compartmental model (see details of the AME derivation in our previous work
[39]). In general, the AME is:

$\displaystyle\frac{d}{dt}x^{\pm}_{k,m,0}=$
$\displaystyle-x^{\pm}_{k,m,0}+\sum_{l}T^{\mp}_{k,m,l}\,x^{\mp}_{k,m,l}$
$\displaystyle-\beta^{\pm}\,(k-m)\,x^{\pm}_{k,m,0}-\gamma^{\pm}\,m\,x^{\mp}_{k,m,0},$
$\displaystyle\frac{d}{dt}x^{\pm}_{k,m,j}=$
$\displaystyle-x^{\pm}_{k,m,j}+A^{\pm}_{k,m,j}\,x^{\pm}_{k,m,j-1}-\beta^{\pm}\,(k-m)\,x^{\pm}_{k,m,j}$
$\displaystyle+\beta^{\pm}\,(k-m+1)\,x^{\pm}_{k,m-1,j-1}$ (3)
$\displaystyle+\gamma^{\pm}\,(m+1)\,x^{\pm}_{k,m+1,j-1}-\gamma^{\pm}\,m\,x^{\pm}_{k,m,j},$

where variables $x^{+}_{k,m,j}$ and $x^{-}_{k,m,j}$ are the fraction of nodes
in state $+1$ or $-1$, respectively, with degree $k$, $m$ neighbors in state
$-1$ that have been $j$ time steps in the current state. The rates
$\beta^{\pm}$ account for the change of state of neighbors ($\pm$) of a node
in state $+1$. The rates $\gamma^{\pm}$ are equivalent but for nodes in state
$-1$. They can be written as

$\displaystyle\beta^{+}=$
$\displaystyle\,\frac{\sum_{j}\sum_{k}p_{k}\sum_{m=0}^{k}(k-m)\,T^{+}_{k,m,j}\,x^{+}_{k,m,j}}{\sum_{j}\sum_{k}p_{k}\sum_{m=0}^{k}(k-m)\,x^{+}_{k,m,j}},$
$\displaystyle\beta^{-}=$
$\displaystyle\,\frac{\sum_{j}\sum_{k}p_{k}\sum_{m=0}^{k}m\,T^{+}_{k,m,j}\,x^{+}_{k,m,j}}{\sum_{j}\sum_{k}p_{k}\sum_{m=0}^{k}m\,x^{+}_{k,m,j}},$
$\displaystyle\gamma^{+}=$
$\displaystyle\,\frac{\sum_{j}\sum_{k}p_{k}\sum_{m=0}^{k}(k-m)\,T^{-}_{k,m,j}\,x^{-}_{k,m,j}}{\sum_{j}\sum_{k}p_{k}\sum_{m=0}^{k}(k-m)\,x^{-}_{k,m,j}},$
(4) $\displaystyle\gamma^{-}=$
$\displaystyle\,\frac{\sum_{j}\sum_{k}p_{k}\sum_{m=0}^{k}m\,T^{-}_{k,m,j}\,x^{-}_{k,m,j}}{\sum_{j}\sum_{k}p_{k}\sum_{m=0}^{k}m\,x^{-}_{k,m,j}},$

where the degree distribution of the chosen network is $p_{k}$. The transition
rate $T^{\pm}_{k,m,j}$ is for the probability of changing state ($\pm\to\mp$)
for an agent of degree $k$, $m$ neighbors in state $-1$ and age $j$, while the
aging rate $A^{\pm}_{k,m,j}$ is for the probability of staying in the same
state and increasing the internal time ($j\to j+1$). For the Symmetrical
Threshold model, these probabilities do not depend on internal time $j$
(Markovian dynamics):

$\displaystyle T^{+}_{k,m,j}=$
$\displaystyle\,\theta\left(\frac{m}{k}-T\right),$ $\displaystyle
T^{-}_{k,m,j}=$ $\displaystyle\,\theta\left(\frac{k-m}{k}-T\right),$ (5)
$\displaystyle A^{\pm}_{k,m,j}=$ $\displaystyle\,1-T^{\pm}_{k,m,j}.$

Figure 3:  Evolution of the average interface density $\rho(t)$ (a), the
average magnetization $m(t)$ (b), the mean internal time $\bar{\tau}(t)$ (c),
and the persistence $p(t)$ (d) for the Symmetrical Threshold model. The
average is computed over $5000$ surviving trajectories (simulations stop when
the system reaches the absorbing ordered states). Results for different values
of $T$ are plotted with diverse markers and colors: red ($T=0.12$) and blue
($T=0.24$) belong to Phase I, green ($T=0.36$) and grey ($T=0.49$) belong to
Phase II and purple ($T=0.86$) belongs to Phase III. Solid colored lines are
the AME integrated solutions, using Eqs. (6)-(8). The initial magnetization is
$m_{0}=0.5$. The system is on an ER graph with $N=4\cdot 10^{4}$ and mean
degree $\langle k\rangle=8$. The dashed green line in (a) shows
$\rho(t)\sim\rho_{0}\,e^{-t}$ and the dashed purple line in (c) shows
$\bar{\tau}(t)=t$ (purple).

If we were not concerned with the internal time dynamics, we can simplify our
AME to the one proposed by J. P. Gleeson in Ref. [48] for general binary-state
models. Here we keep the internal times for a dynamical characterization of
the different phases and as a reference frame for the aging studies in the
next section. The primary approximations of this framework are to assume the
thermodynamic limit ($N\to\infty$) and uncorrelated network with negligible
levels of clustering. For the complex networks considered, these conditions
are satisfied for large $N$, and the differential equations can be solved
numerically using standard methods. The mixed order and ordered frozen
transitions predicted (solid black lines in Figs. 2a and 2b, respectively) are
in agreement with the numerical simulations. The predicted lines represent the
initial and final values of $T$ at which the AME reaches the ordered absorbing
states $m_{f}=\pm 1$. In Fig. 2c, we also observe a good agreement between
numerically integrated solutions (solid colored lines) and numerical
simulations (markers).

An alternative simpler approximation is to consider a heterogeneous mean-field
approximation (HMF) (refer to Appendix A for further details). Although this
approximation captures the qualitative behavior, the numerically integrated
solutions do not agree with numerical simulations (see red dashed lines in
Figs. 2a and 2b, and the colored dotted lines in Fig. 2c), and the frozen
phase is not predicted by this framework. These findings demonstrate that
threshold models need approximations beyond mean-field to achieve accuracy, in
agreement with the findings in Refs. [12, 48, 39].

Beyond the stationary states, the previous phases can be characterized by
their ordering dynamical properties. To describe the coarsening process, we
use the time-dependent average interface density $\rho(t)$ (fraction of links
between nodes in different states), the average magnetization $m(t)$, the mean
internal time $\bar{\tau}(t)$ (mean time spent in the current state over all
the nodes) and the persistence $p(t)$ (fraction of nodes that remain in their
initial state at time $t$) [50]. Fig. 3 shows the average results obtained
from the numerical simulations, starting from an initial magnetization
$m_{0}=0.5$. There are 3 regimes with different dynamical properties:

  * •

Mixed regime (Phase I): It corresponds to Phase I in the static phase diagram
and it is characterized by fast disordering dynamics, which is reflected by an
exponential decay of the persistence. The interface density, the
magnetization, and the mean internal time exhibit fast dynamics towards their
asymptotic values in the dynamically active stationary state (see
$T=0.12,0.24$ in Fig. 3);

  * •

Ordered regime (Phase II): It coincides with Phase II in the static diagram
and it is characterized by an exponential decay of the interface density. The
magnetization tends to the ordered absorbing state based on the initial
majority, and the mean internal time tends to scale as $\bar{\tau}(t)\sim t$.
Persistence in this phase decays until a plateau that corresponds to the
initial majority that reaches consensus (since this fraction of nodes does not
change state from the initial condition). When consensus is reached, the
surviving trajectory is stopped (see $T=0.36,0.49$ in Fig. 3);

  * •

Frozen regime (Phase III): This regime corresponds to Phase III and it is
characterized by an initial ordering process followed by the stop of the
dynamics, with constant values of the metrics. The only exceptions are the
mean internal time that grows as $\bar{\tau}(t)\sim t$ (see $T=0.86$ in Fig.
3) and the persistence.

Using the numerically integrated solutions of AME ($x^{\pm}_{k,m,j}(t)$), we
can compute the magnetization $m(t)$, the interface density $\rho(t)$, and the
mean internal time $\bar{\tau}$:

$\displaystyle\rho(t)=$
$\displaystyle\,\frac{2\,\sum_{j}\sum_{k}p_{k}\sum_{m}mx^{+}_{k,m,j}}{\sum_{j}\sum_{k}p_{k}\sum_{m}k(x^{+}_{k,m,j}+x^{-}_{k,m,j})},$
(6) $\displaystyle m(t)=$
$\displaystyle\,2\sum_{j}\sum_{k}p_{k}\sum_{m}x^{+}_{k,m,j}-1$ (7)
$\displaystyle\bar{\tau}(t)=$
$\displaystyle\,\sum_{j}\sum_{k}p_{k}\sum_{m}j\left(x^{+}_{k,m,j}+x^{-}_{k,m,j}\right).$
(8)

All metrics exhibit a strong agreement between the numerical simulations and
the integrated solutions (see solid lines in Fig. 3). However, the persistence
cannot be directly calculated from the integrated solutions. This is because
the fraction of persistent nodes at time $t$ corresponds to the fraction of
nodes with internal time $j=t$, which is at an extreme of the age distribution
at each time step, since $x^{\pm}_{k,m,j}(t)=0$ for $j>t$. Therefore, the
computation of this measure requires a more sophisticated analysis using
extreme value theory [51].

We note that the dynamical characterization discussed above holds for all
possible $m_{0}$ except for the symmetric initial condition $m_{0}=0$. In this
case, an order-disorder transition arises at a critical mean degree $k_{c}$,
whose value depends on the size of the system $N$ [52].

##  III  Symmetrical Threshold model with aging

Aging refers to the property of agents becoming less likely to change their
state the longer they have remained in that state [30, 44, 36, 37, 35, 34, 38,
39]. In contrast to the original model, which assumes that agents update their
state at a constant rate, this model introduces an activation probability
$p_{A}(j)$ that is inversely proportional to the agent’s internal time $j$. At
each time step, the following two steps are performed:

  1. 1.

A node $i$ with age $j$ is selected at random and activated with probability
$p_{A}(j)$;

  2. 2.

If the fraction of neighbors in a different state is greater than the
threshold $T$, the activated node changes its state from $s_{i}$ to $-s_{i}$
and resets its internal time to $j=0$.

We set the activation probability to $p_{A}(j)=1/(j+2)$ with the aim of
recovering a fat-tailed inter-event time distribution, as observed in simple
contagion models [31, 44].

###  III.1 Mean-field

Figure 4:  Evolution of the average magnetization $m(t)$ (a) and the mean
internal time $\bar{\tau}(t)$ (b) in a complete graph of $N=2500$ nodes.
Results are shown for the Symmetrical Threshold Model (pluses) and the version
with aging (crosses) obtained from simulations. Different colors correspond to
different values of the threshold $T$: red ($T=0.1$) belongs to Phase I, green
($T=0.5$) belongs to Phase II, and purple ($T=0.9$) to Phase III. The initial
magnetization is fixed at $m_{0}=0.5$. The solid and dashed lines correspond
to the numerically integrated solutions from Eq. 11 for the original model
($p_{A}(j)=1$) and the version with aging ($p_{A}(j)=1/(t+2)$), respectively.

Figure 4 compares the evolution of the average magnetization and mean internal
times on a complete graph of the original Symmetrical Threshold model and the
version with aging in phases I, II and III. We observe that, for all
considered threshold values, aging introduces a delay. However, the final
stationary state coincides with the one observed for the original model. To
explain these dynamics, we use a heterogeneous mean-field approach that
considers the effects of aging (HMFA), as in Ref. [34] for other binary-state
models (we use a general HMF description to be applied for a complete graph
and to random networks in next section). In this case, the AME does not work
well due to the high density of the network. For a general network with degree
distribution $p_{k}$, we define the fraction of agents in state $\pm 1$ with
$k$ neighbors and age $j$ at time $t$ as $x^{\pm}_{k,j}(t)$. The probability
of finding a neighbor in state $\pm 1$ is $\tilde{x}^{\pm}$, which can be
written as

$\tilde{x}^{\pm}=\sum_{k}p_{k}\frac{k}{\langle
k\rangle}\,\sum_{j=0}^{\infty}x^{\pm}_{k,j},$ (9)

where $\langle k\rangle$ is the mean degree of the network. The transition
rate $\omega_{k,j}^{\pm}$ for a node with state $\pm 1$, degree $k$ and age
$j$ to change state is given by

$\omega_{k,j}^{\pm}=p_{A}(j)\,\sum_{m=0}^{k}\theta\left(\frac{m}{k}-T\right)\,B_{k,m}[\tilde{x}^{\mp}],$
(10)

where $B_{k,m}[x]$ is the binomial distribution with $k$ attempts, $m$
successes, and with the probability of success $x$. In our model, there are
two possible events for a node with degree $k$ and age $j$:

  * •

It changes state and the age is reset to $j=0$;

  * •

It remains at its state and the age increases by one time step $j=j+1$.

According to these possible events, we can write the rate equations for the
variables $x^{\pm}_{k,j}$ and $x^{\pm}_{k,0}$ as

$\displaystyle\frac{dx^{\pm}_{k,0}}{dt}$
$\displaystyle=\sum_{j=0}^{\infty}x^{\mp}_{k,j}\,\omega_{k,j}^{\mp}-x^{\pm}_{k,0}\qquad\mbox{for
zero age},$ (11) $\displaystyle\frac{dx^{\pm}_{k,j}}{dt}$
$\displaystyle=x^{\pm}_{k,j-1}\,(1-\omega_{k,j-1}^{\pm})-x^{\pm}_{k,j}\qquad\mbox{if
}j>0.$

It can be shown from Eq. (11) that the stationary solution for the fraction of
agents in state $+1$, $x_{f}$, obeys the following implicit equation for a
complete graph (see Appendix B for a detailed explanation):

$x_{f}=\frac{F(1-x_{f})}{F(x_{f})+F(1-x_{f})},$ (12)

where,

$F(x)=1+\sum_{j=1}^{\infty}\prod_{a=0}^{j-1}\left(1-p_{A}(a)\,\sum_{m=(N-1)T}^{N-1}B_{N-1,m}[x]\right).$
(13)

A solution of Eq. (12) can be obtained numerically using standard methods, as
in Ref. [34]. The final magnetization is calculated as $m_{f}=2\,x_{f}-1$.
With this method, we obtain that the phase diagram for the model with aging is
the same as for the original model (refer to Fig. 1a). As a technical point,
we note that a truncation of the summation of the variable $j$ in Eq. (13) is
required for the numerical resolution of the implicit equation. The higher the
maximum age considered $j_{\rm max}$, the higher the accuracy. With $j_{\rm
max}=5\cdot 10^{4}$, the transition lines predicted by this mean-field
approach show great accuracy. Moreover, by numerically integrating Eqs. (11),
the dynamical evolution of the magnetization and mean internal time can be
obtained. Fig. 4 shows the agreement between integrated solutions and Monte
Carlo simulations of the system both for the aging and non-aging versions. It
should be noted that, while aging introduces only a dynamical delay for the
magnetization $m(t)$, the mean internal time $\bar{\tau}(t)$ in Phase I shows
a different dynamical behavior with aging than in the original model. In this
phase, due to the low value of $T$, the agents selected randomly will change
their state (as they fulfill the threshold condition) and reset their internal
time. Consequently, while the internal time fluctuates around a stationary
value for the original model, when aging is incorporated, due to the
activation probability $p_{A}(j)$ chosen, the mean internal time increases
following a recursive relation (Eq. (24)). We refer to Appendix C for a
derivation of this result.

###  III.2  Random networks

Figure 5:  Phase diagram of the Symmetrical Threshold with aging model in an
ER graph of $N=4\cdot 10^{4}$ nodes and $\langle k\rangle=8$. The blue, red,
and green dotted lines show the borders of Phase II (first and last value of
$T$ where the system reaches the absorbing ordered state for each $m_{0}$)
computed from numerical simulations evolving until $t_{\rm max}=10^{3}$,
$10^{4}$ and $10^{5}$ time steps, respectively. Black solid lines show AME
solution integrated $10^{5}$ time steps. Phase ${\rm I}^{*}$, II and III
correspond with the orange, white and gray areas, respectively. The solid
purple lines are the mixed-ordered and ordered-frozen critical lines for the
non-aging version of the model.

In contrast to the results obtained in a complete graph, aging effects have a
significant impact on the phase diagram of the model on random networks. In
Fig. 5, we show the borders of Phase II (first and last value of $T$ where the
system reaches the absorbing ordered state for each $m_{0}$) obtained from
Monte Carlo simulations running up to a maximum time $t_{\rm max}$ (dotted
colored lines). Reaching the stationary state in this model requires a large
number of steps and it has a high computational cost. The two borders of Phase
II exhibit different behavior as we increase the time cutoff $t_{\rm max}$:
while the ordered-frozen border does not change with different $t_{\rm max}$,
the mixed-ordered border is shifted to lower values of $T$ as we increase the
time cutoff $t_{\rm max}$. Our results suggest that Phase I is actually
replaced in a good part of the phase diagram by an ordered phase in which the
absorbing state $m_{f}=\pm 1$ is reached after a large number of time steps.
The ordered-frozen border is now slightly shifted to lower values of the
threshold $T$ due to aging. Similar results are found for a RR graph (see
Appendix D). The dependence of the results with $t_{\rm max}$ calls for a
characterization of different phases in terms of dynamical properties rather
than by the asymptotic value of the magnetization.

Figure 6 shows the time evolution of our ordering metrics. The dynamical
properties are largely affected by the aging mechanism. In terms of the
evolution, we find the following regimes:

  * •

Initial mixing regime (Phase ${\rm{\bf I}}^{*}$): It is characterized by two
dynamical transient regimes: a fast initial disordering dynamics followed by a
slow ordering process. After the initial fast disordering stage, the average
interface density exhibits a very slow (logarithmic-like) decay. Later, due to
the finite size of the system the average interface density follows a power
law decay with time, where $\rho(t)$ scales as $t^{-1}$. This phase exists for
the same domain of parameters ($m_{0}$, $T$) as Phase I (orange region in Fig.
5) in the model without aging (see $T=0.12,0.24$ in Fig. 6);

  * •

Ordered regime (Phase II): It is characterized by a power-law interface decay,
where $\rho(t)$ scales as $t^{-1}$. The magnetization tends to the ordered
absorbing state according to the initial majority (see $T=0.36,0.49$ in Fig.
6);

  * •

Frozen regime (Phase III): It is characterized by an initial tendency towards
the majority consensus, but very fast reaches an absorbing frozen
configuration (see $T=0.86$ in Fig. 6).

Figure 6:  Evolution of the average interface density $\rho(t)$ (a) and the
average magnetization $m(t)$ (b) for the Symmetrical Threshold model with
aging. The average is computed over $5000$ surviving trajectories (simulations
stop when the system reaches the absorbing ordered states) for different
values of $T$, shown by different markers and colors: red ($T=0.12$) and blue
($T=0.24$) belong to Phase ${\rm I}^{*}$, green ($T=0.36$) and grey ($T=0.49$)
belong to Phase II and purple ($T=0.86$) belong to Phase III. The inset in (a)
shows a close look to the evolution for $T=0.12$, in linear-log scale. Solid
colored lines are the AME integrated solutions for $10^{4}$ time steps, using
Eqs. 6 \- 7. The initial magnetization is $m_{0}=0.5$. The system is on an ER
graph with $N=4\cdot 10^{4}$ and mean degree $\langle k\rangle=8$. The dashed
green line in (a) shows $\rho(t)\sim\rho_{0}\,t^{-1}$.

The main effect of aging is that the mixed states of Phase I are no longer
present, at least not for the type of networks that we are analyzing here. We
will show later that Phase I reemerges in denser graphs. Instead, for sparse
graphs, we observe a new Phase ${\rm I}^{*}$ in which the system initially
disorders and later orders until reaching the absorbing states $m_{f}=\pm 1$.
This behavior is shown in Fig. 6 for $T=0.12$ and $0.26$. For $T=0.12$, the
system initially disorders, and then the interface density follows a
logarithmic-like decay (see inset in Fig. 6a). Due to the slow decay, the
system stays in this transient regime even after $10^{6}$ time steps, and the
fall to the absorbing states is not seen. Similarly, for $T=0.26$ the
disordering process stops and then the system gradually evolves towards a
fully ordered state. For this value of $T$, the logarithmic-like decay is not
appreciated and we just observe the power-law decay due to the finite size of
the system. The difference between $T=0.12$ and $T=0.26$ comes from the fact
that in this Phase ${\rm I}^{*}$, the decay of $\rho$ becomes faster as we
increase the threshold $T$ (see inset in Fig. 7). Notice the different
interface decay in Fig. 7 between values of $T<0.3$ (Phase ${\rm I}^{*}$),
where all trajectories show a logarithmic-like decay of $\rho(t)$ in a
transient regime, and $T\geq 0.3$ (Phase II), where trajectories from the
initial condition exhibit ordering dynamics towards the consensus of the
majority. Moreover, we observe that in Phase ${\rm I}^{*}$, the initial
magnetization $m_{0}$ introduces a bias to the stochastic process, implying
that the larger $m_{0}$ in absolute value, the larger the number of
realizations that reach the absorbing state with the same sign of $m_{0}$.
However, the system can still reach the absorbing state of the opposite sign
of $m_{0}$ (initial minority), as shown in the trajectories with $T=0.25$ in
Fig. 7.

Figure 7:  Interface density $\rho(t)$ (lower) and magnetization $m(t)$
(upper) trajectories for different values of the threshold $T$ ($m_{0}=-0.2$)
using the Symmetrical Threshold model with aging. Different colors indicate
different values of $T$. The inset shows a close look at the logarithmic-like
decay, shown in linear-log scale. The system is an ER graph with $N=4\cdot
10^{4}$ and mean degree $\langle k\rangle=8$.

In Phase II, the system asymptotically orders for any initial condition as in
the original model, but the dynamical properties are modified due to the
presence of aging: the exponential decay of the interface density is replaced
by a slow power-law decay, where the exponents of the exponential and the
power-law are found to be the same. Contrary, the dynamical properties of
Phase III are not affected by the presence of aging. The temporal magnitudes
analysis (mean internal time and persistence) can be found in Appendix E.

To account for the results of our Monte Carlo simulations, we use the same
mathematical framework as described in Equation (II.2). According to the
update rules of the Symmetrical Threshold Model with aging, the transition
probabilities now depend on the age $j$, as given by the activation
probability $p_{A}(j)$:

$\displaystyle T^{+}_{k,m,j}=$ $\displaystyle\,p_{A}(j)\,\theta(m/k-T),$
$\displaystyle T^{-}_{k,m,j}=$ $\displaystyle\,p_{A}(j)\,\theta((k-m)/k-T),$
(14) $\displaystyle A^{\pm}_{k,m,j}=$ $\displaystyle\,1-T^{\pm}_{k,m,j}.$

We show in Figure 5 the mixed-ordered and ordered-frozen transition lines
predicted by the integration of the AME equations until a time cutoff $t_{\rm
max}$. We find good agreement between the theoretical predictions and the
simulations both for ER and RR networks (see RR results in Appendix D).
Regarding dynamical properties, the AME integrated solutions exhibit a
remarkable concordance with the evolution of all the metrics as shown in
Figure 6. Minor discrepancies between the numerical simulations and the
integrated solutions can be attributed to the assumption of an infinitely
sized system in the AME.

Figure 8:  Critical threshold $T_{c}$ dependence with the mean degree $\langle
k\rangle$ for the Symmetrical Threshold model with aging for an ER graph with
$N=4\times 10^{4}$ nodes for an initial magnetization of $m_{0}=0.25$ (a) and
$m_{0}=0.75$ (b). The blue and red markers indicate the borders of phases I
and II, which coincide for a sufficiently large value of the mean degree. The
hatched area corresponds to the fulfillment of the inequality in the legend.

The numerical results discussed so far are for random networks with average
degree $\langle k\rangle=8$. According to them and to the analytical insights,
one can conclude that aging significantly changes the phase diagram for sparse
networks. However, we know that the mean-field (fully connected) model with
aging shows the same phase diagram as the model without aging. This implies
that, for ER graphs, as the mean degree $\langle k\rangle$ approaches $N$,
Phase ${\rm I}^{*}$ must disappear. Therefore, the combined effects of
increasing the mean degree and introducing aging need to be investigated in
more detail. Phase II is distinguishable from phases I and ${\rm I}^{*}$
because the system initially orders, i.e., $|\rho_{0}-\rho_{\rm max}|=0$,
where $\rho_{\rm max}$ is the maximum value attained by the interface density
during the dynamical evolution. In contrast, Phase I is distinguished from
Phases ${\rm I}^{*}$ and II because the system remains disordered, i.e.,
$|\rho_{\rm max}-\rho(t_{\rm max})|\approx 0$. Thus, Phase ${\rm I}^{*}$ is
the only phase among these three where $|\rho_{0}-\rho_{\rm max}|>0$ and
$|\rho_{\rm max}-\rho(t_{\rm max})|>0$. Using this criterion, we studied the
dependence of the critical threshold $T_{c}$ on the mean network degree
defining the transition lines between phases I, ${\rm I}^{*}$, and II (see
Fig. 8). In the absence of aging, the red line in Fig. 8 gives the value of
the mixed-ordered transition line $T_{c}$. When aging is included, at low
degree values, Phase I is replaced by ${\rm I}^{*}$, as expected. However, as
the mean degree increases, Phase I emerges despite the presence of aging,
leading to the coexistence of phases I and ${\rm I}^{*}$ in the same phase
diagram over a range of mean degree values. As the mean degree is further
increased, a critical value is reached where Phase ${\rm I}^{*}$ is no longer
present, and the discontinuous transition I-II occurs at the same value than
in the model without aging. Importantly, this critical mean degree at which
Phase ${\rm I}^{*}$ disappears, depends significantly on the initial
magnetization $m_{0}$.

##  IV  Symmetrical Threshold model in a Moore Lattice

We consider next the Symmetrical Threshold model in a Moore lattice, which is
a regular 2-dimensional lattice with interactions among nearest and next-
nearest neighbors ($k=8$). From numerical simulations, we obtain a phase
diagram (Fig. 9a) that is consistent with our previous results in random
networks. The system undergoes a mixed-ordered transition at a threshold value
$T_{c}=2/8$ which is independent of the value of the initial magnetization
$m_{0}$. When $T>4/8$, the system undergoes an ordered-frozen transition at a
critical threshold $T_{c}^{*}$, which depends on $m_{0}$ (similarly to what
happens in random networks). The final magnetization $m_{f}(m_{0})$ (Fig. 9b)
also shows a dependence on $m_{0}$ similar to the one found in RR networks
(Fig. 2c).

Figure 9:  (a) Phase diagram of the Symmetrical Threshold model in a Moore
lattice of size $N=L\times L$, with $L=100$. The color map indicates the value
of the average final magnetization $m_{f}$. Solid black lines are the borders
of Phase II (first and last value of $T$ where the system reaches the
absorbing ordered state for each $m_{0}$), computed from the numerical
simulations. (b) Average final magnetization $m_{f}$ as a function of the
initial magnetization $m_{0}$ for the discrete values of the threshold $T$
(indicated with different colors and markers) in a Moore lattice of the same
size. Average performed over 5000 realizations.

###  IV.1 Original model without aging

Figure 10:  Evolution of the interface density $\rho(t)$ (a), the average
magnetization $m(t)$ (b) and the persistence $p(t)$ (c) for the Symmetrical
Threshold model in a Moore lattice starting from a random configuration with
$m_{0}=0$ (upper) and $m_{0}=0.5$(lower). We plot 50 different trajectories in
solid lines and the average of $5000$ surviving trajectories (simulations stop
when the system reaches the absorbing ordered states) in different markers.
Different colors and markers indicate different threshold values: red
($T=1/8$) and blue ($T=2/8)$ belong to Phase I, green ($T=3/8$) and black
($T=4/8$), and purple ($T=5/8,7/8$) belong to Phase III. The average
magnetization $m(t)$ is computed according to the two symmetric absorbing
states. System size is fixed at $N=L\times L$, $L=200$. The dashed lines in
(a-upper) are $\rho(t)\sim at^{-1/2}$ with $a=0.36$ (black) and $a=0.2$
(green), in (a-lower) are $\rho\sim\exp(-\alpha\cdot t)$ with $\alpha=0.5$
(black) and $\alpha=0.8$ (green), (c) are $p(t)\sim\exp(-\ln(t)^{2})$ (blue).

Fig. 10 shows the results from numerical simulations (for $m_{0}=0$ and $0.5$)
for the average interface density, the magnetization, and the persistence (the
internal time shows the same results as in random graphs). Dynamical
properties change significantly for different values of the threshold and
initial magnetization $m_{0}$. Similarly to the case of random networks, we
find three different regimes corresponding to the three phases, but with some
properties different from the results on random networks:

  * •

Mixed regime (Phase I): It is characterized by fast disordering dynamics with
a persistence decay $p(t)\sim\exp(-\ln(t)^{2})$. The interface density and the
magnetization exhibit fast dynamics towards their asymptotic values in the
dynamically active stationary state (see $T=1/8,2/8$ in Fig. 10);

  * •

Ordered regime (Phase II): It is characterized by an exponential or power-law
decay of the interface density, depending on the initial condition. The
magnetization tends to the absorbing ordered state (see $T=3/8,4/8$ in Fig.
10);

  * •

Frozen regime (Phase III): It is characterized by an initial ordering process,
but the system freezes fast (see $T=5/8$ in Fig. 10).

In particular, in Phase II the persistence and interface density decay as a
power law, $p(t)\sim t^{-0.22}$ and $\rho(t)\sim t^{-1/2}$ for $m_{0}=0$ (as
in Refs. [53, 54, 55, 56]). For a biased initial condition ($m_{0}=0.5$),
$p(t)$ decays to the initial majority fraction (which corresponds to the state
reaching consensus) and $\rho(t)$ follows an exponential decay. Note that, for
$m_{0}=0$, not all trajectories reach the ordered absorbing states ($m_{f}=\pm
1$). There exist other absorbing configurations as, for example, a flat
interface configuration for $T=4/8$, no agent will be able to change, and the
system remains trapped in this state. This result is not observed for
$m_{0}>0$.

Contrary, phases I and III show similar dynamics for a symmetrical ($m_{0}=0$)
and biased ($m_{0}=0.5$) initial conditions. In Phase I, the system shows
disordering dynamics with a persistence decay similar to the one exhibited for
the Voter model in a lattice [50] while in Phase III, the system exhibited
freezing dynamics with an initial tendency towards the majority consensus.

Due to the lattice structure and high clustering, the mathematical tools used
in previous sections for random networks cannot be applied in the case of a
regular lattice. In this case, we restrict ourselves to the results of
numerical simulations.

###  IV.2 The role of aging

Figure 11:  (a) Phase diagram of the Symmetrical Threshold model with aging in
a Moore lattice of $N=L\times L$, with $L=100$. The blue, red and green dotted
lines show the borders of Phase II (first and last value of $T$ where the
system reaches the absorbing ordered state for each $m_{0}$) from numerical
simulations evolving until $t_{\rm max}=10^{3}$, $10^{4}$ and $10^{5}$ time
steps, respectively. Phase ${\rm I}^{*}$, II and III correspond with the
orange, white and gray areas, respectively. The solid purple lines are the
mixed-ordered and ordered-frozen critical lines for the Symmetrical threshold
model (from Fig. 9) (b) Average magnetization at time $t_{\rm max}$
($m_{f}(t_{\rm max})$) as a function of the initial magnetization $m_{0}$ for
different values of the threshold $T$ (indicated with different colors and
markers) in a Moore lattice of $N=L\times L$, with $L=100$. The numerical
simulations are obtained until $t_{\rm max}=10^{4}$ MC steps. Average
performed over 5000 realizations. Figure 12:  Evolution of the average
interface density $\rho(t)$ (a), the average magnetization $m(t)$ (b) and the
persistence $p(t)$ (c) for the Symmetrical model with aging in a Moore lattice
starting from a random configuration with $m_{0}=0$ (upper) and
$m_{0}=0.5$(lower). We plot 30 different trajectories in solid lines and the
average over $5000$ surviving trajectories in symbols. Colors and symbols
indicate different threshold values: red ($T=1/8$) and blue ($T=2/8)$ belong
to Phase ${\rm I}^{*}$, green ($T=3/8$), and black ($T=4/8$) belong to Phase
II, and purple ($T=5/8,7/8$) belong to Phase III. The average magnetization is
computed according to the two symmetric absorbing states. The insets in (a)
show a close look at the evolution for $T=0.12$, in linear-log scale. System
size is fixed at $N=L\times L$, $L=200$. The dashed lines in (a) are $\rho\sim
t^{-\alpha}$ with $\alpha=0.5$ (black) and $\alpha=0.8$ (green), and in (c)
are $p(t)\sim t^{-1}$ (red). Simulations stop when the system reaches the
absorbing ordered states.

We show in Figure 11a borders of Phase II obtained from numerical simulations
running up to a time $t_{\rm max}$ (dotted colored lines). Similarly to the
behavior observed in random networks, the mixed-ordered border is shifted to
lower values of $T$ as we increase the simulation time cutoff $t_{\rm max}$.
Thus, Phase I is replaced by an ordered phase due to the aging mechanism.
Examining the dependence of the final value of the magnetization on its
initial condition $m_{f}(m_{0})$ (Figure 11b), one can conclude that the mixed
phase is still present, at least transiently, as in the initially disordering
phase described in the previous section (Phase ${\rm I}^{*}$). Phase II is
again characterized by an asymptotically ordered state where the initial
majority reaches consensus. However, for this specific structure, near
$m_{0}=0$ and $T=1/2$, the ordered state is not reached for any threshold
value. Furthermore, comparing with Fig. 11b with the results from the model
without aging (Fig. 9b), the discontinuous jump at $m_{0}=0$ for $T=3/8,4/8$
is replaced by a continuous transition, where a range of states with
$0<|m_{f}|<1$ are present around $m_{0}=0$. To determine whether these states
belong to Phase ${\rm I}^{*}$, II or III, we need again a characterization of
phases in terms of dynamical properties. According to the results in Figure
12, we find here the same regimes identified for random networks:

  * •

Initial mixing regime (Phase ${\rm{\bf I}}^{*}$): After the initial
disordering stage, the average interface density shows a very slow decay
reflecting the slow growth of spatial domains in each binary state. The
persistence in this phase shows a power-law decay $p(t)\sim t^{-1}$ (see
$T=1/8,2/8$ in Fig. 12);

  * •

Ordered regime (Phase II): It is characterized by coarsening dynamics that end
in the absorbing states $m_{f}=\pm 1$. The form of the decay of the interface
density depends on the value of $m_{0}$ (see $T=3/8,4/8$ in Fig. 12);

  * •

Frozen regime Phase III): It characterizes by an initial tendency to order but
the system very fast reaches an absorbing frozen configuration (see
$T=5/8,7/8$ in Fig. 12).

The implications of aging become explicit by comparing the dynamical
properties of the cases with aging (Figure 12) and without aging (Figure 10).
When the threshold is $T<3/8$, Phase ${\rm I}$ is replaced by Phase ${\rm
I}^{*}$ in which there is an initial disordering process very fast followed by
a slow coarsening process that accelerates when we increase the threshold.
Although the aging implications in this phase are similar to those observed in
the ER graph, the coarsening process is slower (see insets in Fig. 12a).

In Phase II ($T=3/8,4/8$) and when $m_{0}=0.5$, the system exhibits coarsening
towards the ordered state $m_{f}=\pm 1$. In this case, the exponential decay
$\rho\sim\exp(-\alpha\,t)$ observed in the absence of aging is replaced, due
to aging, by a power law $\rho\sim t^{-\alpha}$ as noted in Ref. [39]. We find
$\alpha=0.5$ and $0.8$ for $T=3/8$ and $4/8$, respectively. For $m_{0}=0$, the
power law decay of the interface density vanishes with aging, and the system
exhibits a coarsening dynamics much slower than for an unbalanced initial
condition. In this region of the phase diagram, spatial clusters start to grow
from the initial condition, but once formed, it takes a long time for the
system to reach the absorbing state $m_{f}=\pm 1$. We note that for these
parameter values, the system is not able to reach $|m|$ over $0.1$ even after
$10^{6}$ time steps, but since there is coarsening from the initial condition,
the expected stationary state as $t\to\infty$ is $m_{f}=\pm 1$. There is
neither initial disordering nor freezing, these values correspond to the
defined Phase II, even though the system exhibits “long-lived segregation” in
a long transient dynamics ( see the difference with the dynamics of the model
without aging in Fig. 13). In Fig. 11a, we differentiate Phase II from Phase
III by analyzing the activity in the system: If agents are changing, even
though the interface decay is slow, the system is in the Phase II. While if
agents are frozen, it lies in Phase III. When comparing the ordered-frozen
critical line to the one from the original model (purple line), we notice that
aging causes certain values ($m_{0}$, $T$) that were previously in Phase II
near the critical line to enter the frozen phase.

Figure 13:  Evolution of a single realization for $T=0.5$ and $m_{0}=0$ using
the Symmetrical threshold model (a) and the version with aging (b). Snapshots
are taken after $1,10,60,440$ and $3300$ time steps in (a) and after
$1,60,3300,2\cdot 10^{5}$ and $5\cdot 10^{6}$ time steps in (b), increasing
from left to right. System size is fixed to $N=L\times L$, $L=256$.

Finally, it should be noted that in Phase ${\rm I}^{*}$, the initial
disordering dynamics drive the system towards $m=0$. Therefore, the subsequent
coarsening dynamics follow the slow interface decay observed in Phase II for
$m_{0}\sim 0$. Thus, the presence of aging implies that the system
asymptotically orders for any initial condition, but due to the initial
disordering, the coarsening dynamics fall into the “long-lived segregation”
regime independently of the initial condition.

##  V  Summary and conclusions

In this work, we have studied with Monte Carlo numerical simulations and
analytical calculations the Symmetrical Threshold Model. In this model, the
agents, nodes of a contact network, can be in one of the two symmetric states
$\pm 1$. System dynamics follows a complex contagion process in which a node
changes state when the fraction of neighboring nodes in the opposite state is
above a given threshold $T$. For $T=1/2$, the model reduces to a majority rule
or the zero temperature Spin Flip Kinetic Ising Model. When the change of
state is only possible in one direction, say from $1$ to $-1$, it reduces to
the Granovetter-Watts Threshold model [6, 9, 39]. We have considered the cases
of a fully connected network, Erdős-Rényi, and random regular networks, as
well as a regular two-dimensional Moore lattice.

We have found that, in the parameter space of threshold $T$ and initial
magnetization $m_{0}$, the model exhibits three distinct phases, namely Phase
${\rm I}$ or mixed, Phase ${\rm II}$ or ordered, and Phase ${\rm III}$ or
frozen. The existence of these three phases is robust for different network
structures. These phases are well characterized by the final state ($m_{f}$),
and by dynamical properties such as the interface density $\rho(t)$, time-
dependent average magnetization $m(t)$, persistence times $p(t)$, and mean
internal time $\bar{\tau}(t)$. These phases can be obtained analytically in
the mean-field case of a fully connected network. For the random networks
considered, we derive an approximate master equation (AME) [48, 39]
considering agents in each state according to their degree $k$, neighbors in
state $-1$, $m$, and age $j$. From this AME, we have also derived a
heterogeneous mean-field (HMF) approximation. While the AME reproduces with
great accuracy the results of Monte Carlo numerical simulations of the model
(both static and dynamic), the HMF shows an important lack of agreement,
highlighting the importance of high-accuracy methods necessary for threshold
models.

Aging is incorporated in the model as a decreasing probability to modify the
state as the time already spent by the agent in that state increases. The key
finding is that the mixed phase (Phase ${\rm I}$), characterized by an
asymptotically disordered dynamically active state, does not always exist: the
aging mechanism can drive the system to an asymptotic absorbing ordered state,
regardless of how low the threshold $T$ is set. A similar effect of aging was
already described for the Schelling model in Ref. [38]. When the dynamics are
examined in detail, a new Phase ${\rm I}^{*}$, defined in terms of dynamical
properties, emerges in the domain of parameters where the model without aging
displays Phase ${\rm I}$. This phase is characterized by an initial
disordering regime ($m\to 0$) followed by a slow ordering dynamics, driving
the system toward the ordered absorbing states (including the one with spins
opposite to the majoritarian initial option). This result is counter-intuitive
since aging incorporates memory into the system, yet in this phase, the system
“forgets” its initial state. The network structure plays an important role in
the emergence of Phase ${\rm I}^{*}$ since it does not exist for complete
graphs. A detailed analysis reveals that Phase ${\rm I}^{*}$ replaces Phase
${\rm I}$ only for sparse networks, including the case of the Moore lattice.
For ER networks we find that, as the mean degree increases, Phase ${\rm I}$
reappears and there is a range of values of the mean degree for which phases
${\rm I}$ and ${\rm I}^{*}$ coexist. Beyond a critical value of the mean
degree, Phase ${\rm I}$ extends over the entire domain of parameters where
Phase ${\rm I}^{*}$ was observed.

While aging favors reaching an asymptotic absorbing ordered state for low
values of $T$ (Phase ${\rm I}$), in Phase II the ordering dynamics are slowed
down by aging, changing, both in random networks and in the Moore lattice, the
exponential decay of the interface density by a power law decay with the same
exponent. The aging mechanism is found not to be important in the frozen Phase
${\rm III}$. All these effects of aging in the three phases are well
reproduced for random networks by the AME derived in this work, which is
general for any chosen activation probability $p_{A}(j)$.

For the Moore lattice, we have also considered in detail the special case of
the initial condition $m_{0}=0$. In this case also Phase ${\rm I}^{*}$ emerges
and Phase ${\rm III}$ is robust against aging effects. However, in Phase ${\rm
II}$ aging destroys the characteristic power law decay of the interface
density $\rho(t)\sim at^{-1/2}$ associated with curvature reduction of domain
walls. This would be a main effect of aging in the dynamics of the phase
transition for the zero temperature spin flip Kinetic Ising model [57].

As a final remark on the general effects of aging in different models of
collective behavior, we note that the replacement of a dynamically active
disordered stationary phase by a dynamically ordering phase is generic. In
this paper, we find the replacement of Phase ${\rm I}$ by Phase ${\rm I}^{*}$.
Likewise in the Voter model, aging destroys long-lived dynamically active
states characterized by a constant value of the average interface density, and
it give rise to coarsening dynamics with a power law decay of the average
interface density [31]. In the same way, in the Schelling segregation model, a
dynamically active mixed phase is replaced, due to the aging effect, by an
ordering phase with segregation in two main clusters. Another aging effect
that seems generic, in phases in which the system orders when there is no
aging, is the replacement of dynamical exponential laws by power laws. This is
what happens here in Phase ${\rm II}$ for the decay of the average interface
density but, likewise, exponential cascades in the Granovetter-Watts model are
replaced due to aging by a power-law growth with the same exponent [39].

Further work with the general AME used in this work would include a new
approach considering the master equation, as in Ref. [58], in order to
incorporate finite size effects (relevant close to $m_{0}=0$) and to give a
mathematical framework to describe the results in Ref. [52].

###### Acknowledgements.

Financial support has been received from the Agencia Estatal de Investigación
(AEI, MCI, Spain) MCIN/AEI/10.13039/501100011033 and Fondo Europeo de
Desarrollo Regional (FEDER, UE) under Project APASOS (PID2021-122256NB-C21 and
PID2021-122256NB-C22) and the María de Maeztu Program for units of Excellence,
grant CEX2021-001164-M.)

##  Appendix A  Heterogeneous mean-field (HMF)

When the transition and aging probabilities do not depend on $j$,
$T^{\pm}_{k,m,j}=T^{\pm}_{k,m}$ and $A^{\pm}_{k,m,j}=A^{\pm}_{k,m}$, if we are
not interested in the solutions $x^{\pm}_{k,m,j}(t)$ and we just want the
final magnetization, Eq. II.2 is reduced to Gleeson’s AME [48] by summing
variable $j$. This is a system of $(k_{\rm max}+1)(k_{\rm max}+1)$
differential equations without loss of accuracy.

Moreover, following the steps in Ref. [48], we perform a heterogeneous mean-
field approximation (HMF) to reduce our system to $k_{\rm max}+1$ differential
equations:

$\displaystyle\frac{d}{dt}x^{-}_{k}=$
$\displaystyle-x^{-}_{k}\sum_{m=0}^{k}T^{-}_{k,m}B_{k,m}[\omega]$
$\displaystyle+\left(1-x^{-}_{k}\right)\sum_{m=0}^{k}T^{+}_{k,m}B_{k,m}[\omega],$
(15)

where $x^{-}_{k}=\sum_{j}\sum_{m}^{k}x^{-}_{k,m,j}$ and
$\omega=\sum_{k}p_{k}\frac{k}{z}x^{-}_{k}$. This system of differential
equations, coupled via $\omega$, cannot be solved analytically. Solving
numerically with standard methods, HMF predicts a mixed-ordered transition
line that qualitatively captures the critical line dependence but
quantitatively differs from the numerical simulations (see the red dashed line
in Figs. 2a and 2b and the dotted colored lines in Fig. 2c). Moreover, this
approximation does not predict a frozen phase in any of the networks
considered. Instead, for high values of $T$, the integrated stationary
solutions are always $m_{f}=\pm 1$, regardless of $m_{0}$. From this analysis,
we conclude that we need sophisticated methods beyond an HMF description to
describe the Symmetrical Threshold model’s phase diagram, as occurs for the
asymmetrical Granovetter-Watts’ Threshold model (see Ref. [39]).

##  Appendix B  Derivation of the stationary solution via the Heterogeneous
mean-field considering the age (HMFA)

Setting the time derivatives to 0 in Eqs. (11), we obtain the relations for
the stationary state:

$\displaystyle x^{\pm}_{k,0}=$
$\displaystyle\sum_{j=0}^{\infty}x^{\mp}_{k,j}\,\omega_{k,j}^{\mp}$
$\displaystyle x^{\pm}_{k,j}=$ $\displaystyle
x^{\pm}_{k,j-1}\,(1-\omega_{k,j-1}^{\pm})\qquad j>0,$ (16)

from where we extract the necessary condition $x^{-}_{k,0}=x^{+}_{k,0}$, as in
Ref. [34]. Notice that by setting $p_{A}(j)=1$ and summing all ages $j$, we
recover the HMF approximation (Eq. 15) for the model without aging. Defining
$x^{\pm}_{j}(t)$ as the fraction of agents in state $\pm 1$ with age $j$:

$x^{\pm}_{j}=\sum_{k}p_{k}\,x^{\pm}_{k,j},$ (17)

and placing the degree distribution of a complete graph $p_{k}=\delta(k-N+1)$
(where $\delta(\cdot)$ is the Dirac delta), we sum variable $k$ and rewrite
Eq. (16) in terms of $x^{\pm}_{j}$:

$\displaystyle x^{\pm}_{0}=$
$\displaystyle\sum_{j=0}^{\infty}x^{\mp}_{j}\,\omega_{j}^{\mp},$
$\displaystyle x^{\pm}_{j}=$ $\displaystyle
x^{\pm}_{j-1}\,(1-\omega_{j-1}^{\pm})\qquad j>0,$ (18)

where $\omega_{j}^{\pm}\equiv\omega_{N-1,j}^{\pm}$. We compute the solution
$x^{\pm}_{j}$ recursively as a function of $x^{\pm}_{0}$:

$x^{\pm}_{j}=x^{\pm}_{0}\,F_{j}^{\pm}\qquad{\rm where}\qquad
F_{j}^{\pm}=\prod_{a=0}^{j-1}(1-\omega_{a}^{\pm}),$ (19)

and summing all $j$,

$x^{\pm}=x^{\pm}_{0}\,F^{\pm}\qquad{\rm where}\qquad
F^{\pm}=1+\sum_{j=1}^{\infty}F_{j}^{\pm}.$ (20)

Using the stationary condition $x^{-}_{0}=x^{+}_{0}$, we reach:

$\frac{x^{+}}{x^{-}}=\frac{F^{+}}{F^{-}}.$ (21)

Notice that, for the complete graph, $\tilde{x}^{+}=x$, $\tilde{x}^{-}=1-x$.
Therefore, $F^{\pm}$ is a function of the variable $x^{\mp}$ ($F^{+}=F(1-x)$).
Thus, we rewrite the previous expression just in terms of the variable $x$:

$\frac{x}{1-x}=\frac{F(1-x)}{F(x)}.$ (22)

##  Appendix C  Internal time recursive relation in Phase ${\rm{\bf
I}}$/${\rm{\bf I}}^{*}$

In Phase I and ${\rm I}^{*}$, the exceeding threshold condition ($m/k>T$) is
full-filled for almost all agents in the system. Thus, agents will change
their state and reset the internal time once activated. For the original
model, all agents are activated once in a time step on average, but for the
model with aging, the activation probability plays an important role. We
consider here a set of $N$ agents that are activated randomly with an
activation probability $p_{A}(j)$ and, once activated, they reset their
internal time. Being $n_{i}(t)$ the fraction of agents with internal time $i$
at the time step $t$, we build a recursive relation for the previously
described dynamics in terms of variables $i$ and $t$:

$\displaystyle n_{1}(t)=\sum_{i=1}^{t-1}p_{A}(i)\,n_{i}(t-1)$ $\displaystyle
n_{i}(t)=(1-p_{A}(i-1))\,n_{i-1}(t-1)\qquad i>1.$ (23)

This recursion relation can be solved numerically from the initial condition
($n_{1}(0)=1$, $n_{i}(0)=0$ for $i>1$). To obtain the mean internal time at
time $t$, we just need to compute the following:

$\bar{\tau}(t)=\sum_{i=1}^{t}i\,n_{i}(t).$ (24)

The solution from this recursive relation describes the mean internal time
dynamics with great agreement with the numerical simulations performed at
Phase I (for the complete graph) and Phase ${\rm I}^{*}$ (for the Erdős-Rényi
and Moore lattice).

##  Appendix D  Symmetrical Threshold model with aging in Random-Regular
graphs

Fig. 14 shows the borders of Phase II (first and last value of $T$ where the
system reaches the absorbing ordered state for each $m_{0}$) obtained from
Monte Carlo simulations running up to a maximum time $t_{\rm max}$ (dotted
colored lines) for a RR graph. Reaching the stationary state in this model
requires a large number of steps and it has a high computational cost. The two
borders of Phase II exhibit different behavior as we increase the maximum
number of time steps $t_{\rm max}$: while the ordered-frozen border does not
change with different $t_{\rm max}$, the mixed-ordered border is shifted to
lower values of $T$ as we increase the simulation time cutoff $t_{\rm max}$.
As it occurs for the results in ER graphs (Fig. 5), our results suggest that
Phase I is actually replaced in a good part of the phase diagram by an ordered
phase in which the absorbing state $m_{f}=\pm 1$ is reached after a large
number of time steps. The ordered-frozen border is now slightly shifted to
lower values of the threshold $T$ due to aging. Figure 14b shows the average
magnetization on RR graphs with simulations running up to a time $t_{\rm
max}=10^{4}$. Upon comparison with Figure 2c, the dependence on $m_{0}$ is
quite similar, indicating the persistence of a transient mixed phase. This
calls for a characterization of different phases in terms of dynamical
properties and not only by the asymptotic value of the magnetization.

Figure 14:  Phase diagram of the Symmetrical Threshold with aging model in a
RR graph (a) of $N=4\cdot 10^{4}$ nodes and $\langle k\rangle=8$. The blue,
red, and green dotted lines show the borders of Phase II (first and last value
of $T$ where the system reaches the absorbing ordered state for each $m_{0}$)
computed from numerical simulations evolving until $t_{\rm max}=10^{3}$,
$10^{4}$ and $10^{5}$ time steps, respectively. Black solid lines show AME
solution integrated $10^{5}$ time steps. Phase ${\rm I}^{*}$, II and III
correspond with the orange, white and gray areas, respectively. The solid
purple lines are the mixed-ordered and ordered-frozen critical lines for the
non-aging version of the model. (b) Average magnetization at time $t_{\rm
max}$ ($m(t_{\rm max})$) as a function of the initial magnetization $m_{0}$
for different values of the threshold $T$ (indicated with different colors and
markers) in an 8-regular graph of $N=4\cdot 10^{4}$. Average performed over
5000 realizations evolved until $t_{\rm max}=10^{4}$ time steps. Dotted and
solid lines are the HMFA (for $T=1/8-4/8$) and AME (for all $T$) solutions
integrated numerically $10^{4}$ time steps.

Regarding to the AME integrated solutions, Figure 14 shows the mixed-ordered
and ordered-frozen transition lines predicted by the integration of the AME
equations until a time cutoff $t_{\rm max}$, which show a good agreement with
the numerical simulations. Figure 14b also shows the predicted dependence of
$m_{f}(m_{0})$ for the RR graph. For comparison purposes, the numerical
integration is computed until the highest $t_{\rm max}$ used in the Monte
Carlo simulations. In addition, we apply the previously introduced HMFA to
these random networks by numerically integrating Eqs. (11). The results,
displayed as dotted colored lines in Figure 14b, show similarity to the AME
solution for $T<0.5$. Nevertheless, as it occurred for the HMF in the original
model, this mathematical framework is not able to describe the frozen phase.

##  Appendix E  Temporal dynamics in the Symmetrical Threshold model with
aging

Figure 15:  Evolution of the mean internal time $\bar{\tau}(t)$ (a) and the
persistence $p(t)$ (b) for the Symmetrical Threshold model with aging. The
average is computed over $5000$ surviving trajectories (simulations stop when
the system reaches the absorbing ordered states) for different values of $T$,
shown by different markers and colors: red ($T=0.12$) and blue ($T=0.24$)
belong to Phase ${\rm I}^{*}$, green ($T=0.36$) and grey ($T=0.49$) belong to
Phase II and purple ($T=0.86$) belong to Phase III. Solid colored lines are
the AME integrated solutions for $10^{4}$ time steps, using Eq. 8. The initial
magnetization is $m_{0}=0.5$. The system is on an Erdős-Rényi graph with
$N=4\cdot 10^{4}$ and mean degree $\langle k\rangle=8$. The dashed lines in
(a) show $\bar{\tau}(t)=t$ (purple) and the solution from the recursive
relation in Eq. (24) (red). The dashed red line in (b) shows $p(t)=t^{-1}$.

Fig. 15 shows the evolution of the temporal dynamics via the mean internal
time and the persistence. The persistence in Phase ${\rm I}^{*}$ shows a
power-law decay, where $p(t)$ scales as $t^{-1}$, and the internal time shows
an increase following the recursive relation given in Equation (24), as it
occurred for the mean-field scenario (Fig. 4). On the other hand, in Phase II,
the persistence decays from $1$ to the fraction of nodes of the initial
majority (the one that does not change state and reaches consensus) and the
mean internal time scales linearly with time, $\bar{\tau}(t)\sim t$. For the
internal time, the AME integrated solutions exhibit a remarkable concordance
with the numerical simulations. Minor discrepancies between the numerical
simulations and the integrated solutions can be attributed to the assumption
of an infinitely sized system in the AME. As it occurred for the model without
aging, the persistence cannot be predicted by this framework.

## References

  * Liggett  _et al._ [1999] T. M. Liggett _et al._ ,  _Stochastic interacting systems: contact, Voter and exclusion processes_ , Vol. 324 (springer science & Business Media, 1999).
  * Sood and Redner [2005] V. Sood and S. Redner, Voter Model on Heterogeneous Graphs, Physical Review Letters 94, 10.1103/physrevlett.94.178701 (2005).
  * Suchecki  _et al._ [2005] K. Suchecki, V. M. Eguíluz, and M. San Miguel, Voter model dynamics in complex networks: Role of dimensionality, disorder, and degree distribution, Phys. Rev. E 72, 036132 (2005).
  * Fernández-Gracia  _et al._ [2014] J. Fernández-Gracia, K. Suchecki, J. J. Ramasco, M. San Miguel, and V. M. Eguíluz, Is the Voter Model a Model for Voters?, Physical Review Letters 112, 10.1103/physrevlett.112.158701 (2014).
  * Redner [2019] S. Redner, Reality-inspired Voter models: A mini-review, Comptes Rendus Physique 20, 275 (2019).
  * Granovetter [1978] M. Granovetter, Threshold Models of Collective Behavior, American Journal of Sociology 83, 1420 (1978).
  * Pastor-Satorras  _et al._ [2015] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Epidemic processes in complex networks, Reviews of Modern Physics 87, 925 (2015).
  * Gleeson [2011] J. P. Gleeson, High-Accuracy Approximation of Binary-State Dynamics on Networks, Physical Review Letters 107, 10.1103/physrevlett.107.068701 (2011).
  * Watts [2002] D. J. Watts, A simple model of global cascades on random networks, Proceedings of the National Academy of Sciences 99, 5766 (2002).
  * Centola  _et al._ [2007] D. Centola, V. M. Eguíluz, and M. W. Macy, Cascade dynamics of complex propagation, Physica A: Statistical Mechanics and its Applications 374, 449 (2007).
  * unk [2018] How Behavior Spreads: The Science of Complex Contagions How Behavior Spreads: The Science of Complex Contagions Damon Centola Princeton University Press, 2018. 308 pp., Science 361, 1320 (2018).
  * Gleeson and Cahalane [2007] J. P. Gleeson and D. J. Cahalane, Seed size strongly affects cascades on random networks, Physical Review E 75, 10.1103/physreve.75.056103 (2007).
  * Hackett  _et al._ [2011] A. Hackett, S. Melnik, and J. P. Gleeson, Cascades on a class of clustered random networks, Physical Review E 83, 10.1103/physreve.83.056107 (2011).
  * Hackett and Gleeson [2013] A. Hackett and J. P. Gleeson, Cascades on clique-based graphs, Physical Review E 87, 10.1103/physreve.87.062801 (2013).
  * Gleeson [2008] J. P. Gleeson, Cascades on correlated and modular random networks, Physical Review E 77, 10.1103/physreve.77.046117 (2008).
  * de Arruda  _et al._ [2020] G. F. de Arruda, G. Petri, and Y. Moreno, Social contagion models on hypergraphs, Physical Review Research 2, 10.1103/physrevresearch.2.023032 (2020).
  * Diaz-Diaz  _et al._ [2022] F. Diaz-Diaz, M. San Miguel, and S. Meloni, Echo chambers and information transmission biases in homophilic and heterophilic networks, Scientific Reports 12, 10.1038/s41598-022-13343-6 (2022).
  * Min and San Miguel [2023] B. Min and M. San Miguel, Threshold cascade dynamics in coevolving networks, Entropy 25, 929 (2023).
  * de Oliveira [1992] M. J. de Oliveira, Isotropic majority-vote model on a square lattice, Journal of Statistical Physics 66, 273 (1992).
  * Pereira and Moreira [2005] L. F. Pereira and F. B. Moreira, Majority-vote model on random graphs, Physical Review E 71, 016123 (2005).
  * Campos  _et al._ [2003] P. R. Campos, V. M. de Oliveira, and F. B. Moreira, Small-world effects in the majority-vote model, Physical Review E 67, 026104 (2003).
  * Castellano  _et al._ [2009a] C. Castellano, M. A. Muñoz, and R. Pastor-Satorras, Nonlinearq-Voter model, Physical Review E 80, 10.1103/physreve.80.041129 (2009a).
  * Mobilia [2015] M. Mobilia, Nonlinear q-voter model with inflexible zealots, Physical Review E 92, 012803 (2015).
  * Mellor  _et al._ [2016] A. Mellor, M. Mobilia, and R. Zia, Characterization of the nonequilibrium steady state of a heterogeneous nonlinear q-voter model with zealotry, EPL (Europhysics Letters) 113, 48001 (2016).
  * Min and San Miguel [2017] B. Min and M. San Miguel, Fragmentation transitions in a coevolving nonlinear voter model, Scientific Reports 7, 12864 (2017).
  * Jedrzejewski [2017] A. Jedrzejewski, Pair approximation for the q-voter model with independence on complex networks, Physical Review E 95, 012307 (2017).
  * Peralta  _et al._ [2018] A. F. Peralta, A. Carro, M. San Miguel, and R. Toral, Analytical and numerical study of the non-linear noisy Voter model on complex networks, Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 075516 (2018).
  * Nowak and Sznajd-Weron [2019] B. Nowak and K. Sznajd-Weron, Homogeneous symmetrical threshold model with nonconformity: Independence versus anticonformity, Complexity 2019 (2019).
  * Nowak and Sznajd-Weron [2020] B. Nowak and K. Sznajd-Weron, Symmetrical threshold model with independence on random graphs, Physical Review E 101, 052316 (2020).
  * Stark  _et al._ [2008] H.-U. Stark, C. J. Tessone, and F. Schweitzer, Decelerating Microdynamics Can Accelerate Macrodynamics in the Voter Model, Physical Review Letters 101, 10.1103/physrevlett.101.018701 (2008).
  * Fernández-Gracia  _et al._ [2011] J. Fernández-Gracia, V. M. Eguíluz, and M. San Miguel, Update rules and interevent time distributions: Slow ordering versus no ordering in the Voter model, Physical Review E 84, 10.1103/physreve.84.015103 (2011).
  * Pérez  _et al._ [2016] T. Pérez, K. Klemm, and V. M. Eguíluz, Competition in the presence of aging: dominance, coexistence, and alternation between states, Scientific Reports 6, 10.1038/srep21128 (2016).
  * Boguñá  _et al._ [2014] M. Boguñá, L. F. Lafuerza, R. Toral, and M. A. Serrano, Simulating non-markovian stochastic processes, Phys. Rev. E 90, 042108 (2014).
  * Chen  _et al._ [2020] H. Chen, S. Wang, C. Shen, H. Zhang, and G. Bianconi, Non-Markovian majority-vote model, Physical Review E 102, 10.1103/physreve.102.062311 (2020).
  * Peralta  _et al._ [2020a] A. F. Peralta, N. Khalil, and R. Toral, Ordering dynamics in the Voter model with aging, Physica A: Statistical Mechanics and its Applications 552, 122475 (2020a).
  * Artime  _et al._ [2018] O. Artime, A. F. Peralta, R. Toral, J. J. Ramasco, and M. San Miguel, Aging-induced continuous phase transition, Physical Review E 98, 10.1103/physreve.98.032104 (2018).
  * Peralta  _et al._ [2020b] A. F. Peralta, N. Khalil, and R. Toral, Reduction from non-Markovian to Markovian dynamics: the case of aging in the noisy-Voter model, Journal of Statistical Mechanics: Theory and Experiment 2020, 024004 (2020b).
  * Abella  _et al._ [2022] D. Abella, M. San Miguel, and J. J. Ramasco, Aging effects in schelling segregation model, Scientific Reports 12, 10.1038/s41598-022-23224-7 (2022).
  * Abella  _et al._ [2023] D. Abella, M. San Miguel, and J. J. Ramasco, Aging in binary-state models: The threshold model for complex contagion, Physical Review E 107, 024101 (2023).
  * Iribarren and Moro [2009] J. L. Iribarren and E. Moro, Impact of Human Activity Patterns on the Dynamics of Information Diffusion, Physical Review Letters 103, 10.1103/physrevlett.103.038702 (2009).
  * Karsai  _et al._ [2011] M. Karsai, M. Kivelä, R. K. Pan, K. Kaski, J. Kertész, A.-L. Barabási, and J. Saramäki, Small but slow world: How network topology and burstiness slow down spreading, Physical Review E 83, 10.1103/physreve.83.025102 (2011).
  * Rybski  _et al._ [2012] D. Rybski, S. V. Buldyrev, S. Havlin, F. Liljeros, and H. A. Makse, Communication activity in a social network: relation between long-term correlations and inter-event clustering, Scientific Reports 2, 10.1038/srep00560 (2012).
  * Zignani  _et al._ [2016] M. Zignani, A. Esfandyari, S. Gaito, and G. P. Rossi, Walls-in-one: usage and temporal patterns in a social media aggregator, Applied Network Science 1, 10.1007/s41109-016-0009-9 (2016).
  * Artime  _et al._ [2017] O. Artime, J. J. Ramasco, and M. San Miguel, Dynamics on networks: competition of temporal and topological correlations, Scientific Reports 7, 10.1038/srep41627 (2017).
  * Kumar  _et al._ [2020] P. Kumar, E. Korkolis, R. Benzi, D. Denisov, A. Niemeijer, P. Schall, F. Toschi, and J. Trampert, On interevent time distributions of avalanche dynamics, Scientific Reports 10, 10.1038/s41598-019-56764-6 (2020).
  * Erdős  _et al._ [1960] P. Erdős, A. Rényi, _et al._ , On the evolution of random graphs, Publ. Math. Inst. Hung. Acad. Sci 5, 17 (1960).
  * Wormald  _et al._ [1999] N. C. Wormald _et al._ , Models of random regular graphs, London Mathematical Society Lecture Note Series , 239 (1999).
  * Gleeson [2013] J. P. Gleeson, Binary-State Dynamics on Complex Networks: Pair Approximation and Beyond, Physical Review X 3, 10.1103/physrevx.3.021004 (2013).
  * Castellano  _et al._ [2009b] C. Castellano, S. Fortunato, and V. Loreto, Statistical physics of social dynamics, Rev. Mod. Phys. 81, 591 (2009b).
  * Ben-Naim  _et al._ [1996] E. Ben-Naim, L. Frachebourg, and P. L. Krapivsky, Coarsening and persistence in the voter model, Physical Review E 53, 3078 (1996).
  * Haan and Ferreira [2006] L. Haan and A. Ferreira,  _Extreme value theory: an introduction_ , Vol. 3 (Springer, 2006).
  * Pournaki  _et al._ [2022] A. Pournaki, E. Olbrich, S. Banisch, and K. Klemm, Order-disorder transition in the zero-temperature ising model on random graphs (2022).
  * Stauffer  _et al._ [1994] D. Stauffer, J. Adler, and A. Aharony, Universality at the three-dimensional percolation threshold, Journal of Physics A: Mathematical and General 27, L475 (1994).
  * Derrida [1995] B. Derrida, Exponents appearing in the zero-temperature dynamics of the 1D Potts model, Journal of Physics A: Mathematical and General 28, 1481 (1995).
  * Derrida  _et al._ [1995] B. Derrida, V. Hakim, and V. Pasquier, Exact First-Passage Exponents of 1D Domain Growth: Relation to a Reaction-Diffusion Model, Physical Review Letters 75, 751 (1995).
  * Derrida [1997] B. Derrida, How to extract information from simulations of coarsening at finite temperature, Physical Review E 55, 3705 (1997).
  * Gunton [1983] J. Gunton, The dynamics of first-order phase transitions, Phase transitions and critical phenomena 8, 267 (1983).
  * Peralta and Toral [2020] A. F. Peralta and R. Toral, Binary-state dynamics on complex networks: Stochastic pair approximation and beyond, Physical Review Research 2, 10.1103/physrevresearch.2.043370 (2020).