Ordering dynamics and aging in the Symmetrical Threshold model

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.


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 binarystate 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 twodimensional 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) [39,48], 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
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 = ±1.Every agent has assigned a fixed threshold 0 ≤ T ≤ 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 → −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 ) i s i , fluctuates around a constant value.

A. 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 = ±1) according to the initial magnetization m 0 ; • Phase III or Frozen: The system freezes at the initial random state m f = m 0 .
For a given initial magnetization m 0 = 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 θ(x − T ), where θ is the Heaviside step function.Thus, in the thermodynamic limit (N → ∞), the variable x evolves over time according to the following mean-field equation: Here, V (x) is the potential function.The stationary value of x, x st , is the solution of the implicit equation resulting from setting the time derivative equal to 0. No analytical expression of x st has been found in terms of T , but the solutions can be understood in terms of the potential V (x): 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. ( 2).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 = ±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 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 [1,3,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 = ±1.

B. 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 k = 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 mag-netization 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 ±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: 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 β ± account for the change of state of neighbors (±) of a node in state +1.The rates γ ± are equivalent but for nodes in state −1.They can be written as , where the degree distribution of the chosen network is p k .The transition rate T ± k,m,j is for the probability of changing state (± → ∓) for an agent of degree k, m neighbors in state −1 and age j, while the aging rate A ± k,m,j is for the probability of staying in the same state and increasing the internal time (j → j + 1).For the Symmetrical Threshold model, these probabilities do not depend on internal time j (Markovian dynamics): 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 → ∞) 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 = ±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,39,48].
Beyond the stationary states, the previous phases can be characterized by their ordering dynamical properties.To describe the coarsening process, we use the timedependent average interface density ρ(t) (fraction of links between nodes in different states), the average magnetization m(t), the mean internal time τ (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 τ (t) ∼ 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 τ (t) ∼ t (see T = 0.86 in Fig. 3) and the persistence.
Using the numerically integrated solutions of AME (x ± k,m,j (t)), we can compute the magnetization m(t), the interface density ρ(t), and the mean internal time τ : 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 ± 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,[34][35][36][37][38][39]44].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.A node i with age j is selected at random and activated with probability p A (j); 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].

A. Mean-field
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 ±1 with k neighbors and age j at time t as x ± k,j (t).The probability of finding a neighbor in state ±1 is x± , which can be written as where k is the mean degree of the network.The transition rate ω ± k,j for a node with state ±1, degree k and age j to change state is given by 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 ± k,j and x ± k,0 as for zero age, (11) 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): where, (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 max , the higher the accuracy.With j max = 5 • 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 τ (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.(C2)).We refer to Appendix C for a derivation of this result.

B. Random networks
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 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 max : while the orderedfrozen border does not change with different t max , the 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 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 ρ(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 ρ(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).
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 I * in which the system initially disorders and later orders until reaching the absorbing states m f = ±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 I * , the decay of ρ 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 I * ), where all trajectories show a logarithmic-like decay of ρ(t) in a transient regime, and T ≥ 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 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.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 (3).According to the update rules of the Symmetrical Threshold Model with aging, the tran-sition probabilities now depend on the age j, as given by the activation probability p A (j): A ± k,m,j = 1 − T ± k,m,j .We show in Figure 5 the mixed-ordered and orderedfrozen transition lines predicted by the integration of the AME equations until a time cutoff t 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.
The numerical results discussed so far are for random networks with average degree k = 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 k approaches N , Phase 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 I * because the system initially orders, i.e., |ρ 0 − ρ max | = 0, where ρ max is the maximum value attained by the interface density during the dynamical evolution.In contrast, Phase I is distinguished from Phases I * and II because the system remains disordered, i.e., work degree defining the transition lines between phases I, 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 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 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 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 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). A. Original model without aging 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:

b a
• Mixed regime (Phase I): It is characterized by fast disordering dynamics with a persistence decay p(t) ∼ 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) ∼ t −0.22 and ρ(t) ∼ 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 ρ(t) follows an exponential decay.Note that, for m 0 = 0, not all trajectories reach the ordered absorbing states (m f = ±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.

B. The role of aging
We show in Figure 11a   (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 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 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 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 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) ∼ 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 = ±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 (Fig- ure 12) and without aging (Figure 10).When the threshold is T < 3/8, Phase I is replaced by Phase 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 = ±1.In this case, the exponential decay ρ ∼ exp(−α t) observed in the absence of aging is replaced, due to aging, by a power law ρ ∼ t −α as noted in Ref. [39].We find α = 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 = ±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 → ∞ is m f = ±1.There is neither initial disordering nor freezing, these values correspond to the defined Phase II, even though the system exhibits "longlived 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 orderedfrozen 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.
Finally, it should be noted that in Phase 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 ∼ 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 ±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 I or mixed, Phase II or ordered, and Phase 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 ρ(t), time-dependent average magnetization m(t), persistence times p(t), and mean internal time τ (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) [39,48] 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 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 I * , defined in terms of dynamical properties, emerges in the domain of parameters where the model without aging displays Phase I.This phase is characterized by an initial disordering regime (m → 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 I * since it does not exist for complete graphs.A detailed analysis reveals that Phase I * replaces Phase I only for sparse networks, including the case of the Moore lattice.For ER networks we find that, as the mean degree increases, Phase I reappears and there is a range of values of the mean degree for which phases I and I * coexist.Beyond a critical value of the mean degree, Phase I extends over the entire domain of parameters where Phase I * was observed.
While aging favors reaching an asymptotic absorbing ordered state for low values of T (Phase 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 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 I * emerges and Phase III is robust against aging effects.However, in Phase II aging destroys the characteristic power law decay of the interface density ρ(t) ∼ 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 I by Phase 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 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].and placing the degree distribution of a complete graph p k = δ(k − N + 1) (where δ(•) is the Dirac delta), we sum variable k and rewrite Eq. (B1) in terms of x ± j : x where ω ± j ≡ ω ± N −1,j .We compute the solution x ± j recursively as a function of x ± 0 : x and summing all j, Using the stationary condition x − 0 = x + 0 , we reach: Notice that, for the complete graph, x+ = x, x− = 1 − x.Therefore, F ± is a function of the variable x ∓ (F + = F (1 − x)).Thus, we rewrite the previous expression just in terms of the variable x: In Phase I and 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: This recursion relation can be solved numerically from 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: 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 I * (for the Erdős-Rényi and Moore lattice).maximum number of time steps t max : while the orderedfrozen border does not change with different t max , the mixed-ordered border is shifted to lower values of T as we increase the simulation time cutoff t 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 = ±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 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.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 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 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 HMF in the original model, this mathematical framework is not able to describe the frozen phase.Fig. 15 shows the evolution of the temporal dynamics via the mean internal time and the persistence.The persistence in Phase 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 (C2), 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, τ (t) ∼ 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.
FIG. 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 − |m0|)/2 and T = (1 + |m0|)/2, respectively.Average performed over 5000 realizations.(b) Potential representation from Eq. (2) for a set of values of the threshold T , shown in different colors.

FIG. 2 .
FIG. 2. Phase diagram of the Symmetrical Threshold model in an ER (a) and a RR (b) graph, both of N = 4 • 10 4 nodes and mean degree k = 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 m0 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.

FIG. 3 .
FIG. 3. Evolution of the average interface density ρ(t) (a), the average magnetization m(t) (b), the mean internal time τ (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 m0 = 0.5.The system is on an ER graph with N = 4 • 10 4 and mean degree k = 8.The dashed green line in (a) shows ρ(t) ∼ ρ0 e −t and the dashed purple line in (c) shows τ (t) = t (purple).

FIG. 4 .
FIG. 4. Evolution of the average magnetization m(t) (a) and the mean internal time τ (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 m0 = 0.5.The solid and dashed lines correspond to the numerically integrated solutions from Eq. 11 for the original model (pA(j) = 1) and the version with aging (pA(j) = 1/(t + 2)), respectively.

FIG. 5 .
FIG. 5. Phase diagram of the Symmetrical Threshold with aging model in an ER graph of N = 4 • 10 4 nodes and k = 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 m0) computed from numerical simulations evolving until tmax = 10 3 , 10 4 and 10 5 time steps, respectively.Black solid lines show AME solution integrated 10 5 time steps.Phase 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.

FIG. 6 .
FIG. 6. Evolution of the average interface density ρ(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 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 m0 = 0.5.The system is on an ER graph with N = 4 • 10 4 and mean degree k = 8.The dashed green line in (a) shows ρ(t) ∼ ρ0 t −1 .

FIG. 7 .
FIG. 7. Interface density ρ(t) (lower) and magnetization m(t) (upper) trajectories for different values of the threshold T (m0 = −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•10 4 and mean degree k = 8.

FIG. 8 .
FIG. 8. Critical threshold Tc dependence with the mean degree k for the Symmetrical Threshold model with aging for an ER graph with N = 4 × 10 4 nodes for an initial magnetization of m0 = 0.25 (a) and m0 = 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.

FIG. 9 .
FIG. 9. (a) Phase diagram of the Symmetrical Threshold model in a Moore lattice of size N = L × 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 m0), computed from the numerical simulations.(b) Average final magnetization m f as a function of the initial magnetization m0 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.
FIG. 11.(a) Phase diagram of the Symmetrical Threshold model with aging in a Moore lattice of N = L × 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 m0) from numerical simulations evolving until tmax = 10 3 , 10 4 and 10 5 time steps, respectively.Phase 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 tmax (m f (tmax)) as a function of the initial magnetization m0 for different values of the threshold T (indicated with different colors and markers) in a Moore lattice of N = L × L, with L = 100.The numerical simulations are obtained until tmax = 10 4 MC steps.Average performed over 5000 realizations.

FIG. 12 .
FIG. 12. Evolution of the average interface density ρ(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 m0 = 0 (upper) and m0 = 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 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 × L, L = 200.The dashed lines in (a) are ρ ∼ t −α with α = 0.5 (black) and α = 0.8 (green), and in (c) are p(t) ∼ t −1 (red).Simulations stop when the system reaches the absorbing ordered states.

FIG. 13 .
FIG. 13.Evolution of a single realization for T = 0.5 and m0 = 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 • 10 5 and 5 • 10 6 time steps in (b), increasing from left to right.System size is fixed to N = L × L, L = 256.

)
Appendix C: Internal time recursive relation in Phase I/I *

FIG. 14 .
Fig.14shows 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 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

FIG. 15 .
FIG. 15.Evolution of the mean internal time τ (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 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 m0 = 0.5.The system is on an Erdős-Rényi graph with N = 4 • 10 4 and mean degree k = 8.The dashed lines in (a) show τ (t) = t (purple) and the solution from the recursive relation in Eq. (C2) (red).The dashed red line in (b) shows p(t) = t −1 .
Appendix E: Temporal dynamics in the SymmetricalThreshold model with aging