Influence of cumulative damage on synchronization of Kuramoto oscillators on networks

In this paper, we study the synchronization of identical Kuramoto phase oscillators under cumulative stochastic damage to the edges of networks. We analyze the capacity of coupled oscillators to reach a coherent state from initial random phases. The process of synchronization is a global function performed by a system that gradually changes when the damage weakens individual connections of the network. We explore diverse structures characterized by different topologies. Among these are deterministic networks as a wheel or the lattice formed by the movements of the knight on a chess board, and random networks generated with the Erd\H{o}s-R\'enyi and Barab\'asi-Albert algorithms. In addition, we study the synchronization times of 109 non-isomorphic graphs with six nodes. The synchronization times and other introduced quantities are sensitive to the impact of damage, allowing us to measure the reduction of the capacity of synchronization and classify the effect of damage in the systems under study. This approach is general and paves the way for the exploration of the effect of damage accumulation in diverse dynamical processes in complex systems.


Introduction
Synchronization is a collective process in which a coupled population, under certain conditions, becomes self-organized in such a way that their components evolve to follow the same dynamical pattern [1,2,3,4].This process is one of the most attractive phenomena in nature and some common examples are the synchronization of flashing fireflies [2,5], a crowd clapping in unison [6], synchronization in arrays of Josephson junctions [7], among many others [2,3,8].Synchronization process is considered universal [2] so it is present in extremely diverse systems and plays a fundamental role in their functioning.In particular, synchronous activity is pivotal in groups of living beings such as colonies of insects, the flock of birds, or the school of fish, since by coherent behavior they can coordinate their movement to find food or to face threats [9].Also, it is involved in vital functions e.g circadian rhythms [10,11], the functioning of the heart due to the synchronization of pacemaker cells [12,13] and, physiological brain activity [14,15,16].Conversely, abnormal neural synchronization is linked to neurological disorders like epilepsy or Parkinson's disease [17,18].
On the other hand, diverse complex systems suffer a gradual reduction in their capacity to perform specific functions due to the accumulation of damage.Although most of them possess repairing mechanisms, when the damage is severe this reparation process can be incomplete or incorrect due to the urgency of the system to recover without compromising its global performance [42].A good example of this is the appearance of scars after an injury [43].The remaining damage or "misrepairs" accumulates over time and degrades the functionality of the system [42].In the case of complex structures, such as living beings, this phenomenon of damage accumulation can be understood as the mechanism that generates aging.One theory proposed to explain aging in living beings as the accumulation of residual damage resulting from imperfect repair processes, spanning from molecules to tissues.As damage accumulates, it progressively impairs vital functions, making them more vulnerable to environmental hazards, increasing the potential risk of disease, and eventually culminating in death [42].The whole process of accumulation of damage and degradation of systems can be observed in social systems and civilizations [44], or in the reduction of the transport capacity in a complex system [43,45,46], just to mention a few examples.
The study of dynamical processes in systems with damage accumulation to establish global measures that quantify the degradation of their functionality is important to identify vulnerabilities in systems or to assess the performance of a process.Also, it helps to better understand the relation between the dynamics and the structure of a system to develop new designs robust and resilient to damage.As we mentioned, synchronization process is a universal phenomenon; in particular, its role in vital functions makes it attractive to study in the context of a system with damage accumulation.For example, the synapses between neurons are responsible for the connections in neural circuits, but some dysfunctions can alter their efficiency [47].Likewise, intracellular communication may be restricted by aging, potentially influencing the synchronization process of pacemaker cells in the heart [48].Across various systems, such as power grids, a similar situation arises where transmission lines frequently encounter reduced capacity.This outcome can be attributed to factors like inadequate maintenance, unfavorable weather conditions, and contamination, among others [49].Some researchers have investigated synchronization process in systems that are exposed to severe damage, which implies removing some components of the system.In network science, this effect is modeled with the complete removal of particular groups of links or nodes [26,50,51,52,53,54,55], where the percolation connectivity threshold defines the limit of the global functionality of the system.Furthermore, different authors have explored synchronization in timevarying topologies with modifications in the network structure occurring at a temporal scale comparable to the characteristic time scale of synchronization [56,57,58,59,60].A recent work included perturbations in the variables used to characterize the dynamics of oscillators [61].As far as we know, synchronization processes in systems that evolve with gradual accumulation of damage have been explored less.
In this contribution, we study the influence of cumulative stochastic damage in the synchronization of Kuramoto oscillators on networks.To this end, we use identical phase Kuramoto oscillators that are coupled through a weighted network where links represent the local coupling of the oscillators, and the accumulation of damage is modeled with a preferential attachment mechanism where edges with previous faults are more likely to be damaged.Also, the damage affects the structure in a non-symmetric way and decreases the coupling strength of the oscillators, provoking alterations in the capacity of the system to reach synchronization from random phases.We analyze different structures that fully synchronize in the absence of any damage and have different topologies.We include two deterministic networks as a wheel and the lattice formed by the movements of the knight on a chess board that contrast with two random networks generated with the Erdős-Rényi and Barabási-Albert algorithms.Additionally, we analyze 109 non-isomorphic networks with six nodes.Our findings allow us to classify the effect of damage in networks.The explored framework is general and introduces different tools for the characterization of the effect of accumulation of damage in diverse dynamical processes in complex systems.

Preliminaries
We consider a system with N nodes formed by identical Kuramoto oscillators with a coupling structure defined by a network where the natural frequency ω is the same for each oscillator [62].However, all the methods developed can be adapted to the analysis of systems where the oscillators have different characteristics; for example, in systems where natural frequencies are sampled from a probability distribution.At a time t, the oscillators are characterized by their phases θ i (t) with i = 1, 2, . . ., N [21,22,40].The oscillators in the original Kuramoto model are all-to-all coupled [19]; however, the coupling can be defined using the connectivity information of a network [21,22].This coupling network corresponds to a simple connected graph formed by a set of N nodes denoted by V and a set of edges E that connect pairs of nodes (i, j), |E| is the total number of different edges in the network.The topology of the structure is described by its N × N adjacency matrix A with elements A ij = A ji = 1 if the nodes i, j are connected and A ij = 0 otherwise; the diagonal elements are A ii = 0 since the nodes are not connected with themselves.In this system, the evolution of phases of identical Kuramoto oscillators placed in the nodes is described by the system of nonlinear coupled differential equations [62] for i = 1, 2, . . ., N and, where K > 0 corresponds to the global coupling strength of the system.In addition, the system in Eq. ( 1) satisfies rotational symmetry so it is invariant under the transformation θ i → θ i + ωt.In this way, it is possible to set ω = 0 and rescale the time by setting K = 1 [63].Then, Eq. ( 1) is transformed into [62,63,64] The Kuramoto system in Eq. ( 2) is a gradient system that eventually evolves to a fixed point [62,63,64].Some of the equilibrium states correspond to coherent states, which means all the oscillators in the system take the same phase value.In other words, the system is completely phase synchronized [62].The conditions under which a Kuramoto system evolves to complete synchronization are still under study but the topology of the system plays an important role [62,63,65,66].A way to assess the phase coherence of the Kuramoto oscillators is through the macroscopic order parameter [21,22,40] where i = √ −1.From Eq. (3), 0 ≤ r(t) ≤ 1.In particular, if the Kuramoto system exhibits complete phase coherence r(t) = 1, on the contrary r(t) = 0 if the oscillators  2) using random initial phases sampled from a uniform distribution in the interval [0, 2π), the value of the phases is codified in the color bar.(c) Evolution of the order parameter r(t) in the synchronization of the network.Thin lines represent the results of 1000 realizations using random initial phases and the dashed line shows the ensemble average denoted as r(t) .move in a completely incoherent manner or form clusters of synchronized oscillators.
As an example of a Kuramoto system that evolves to a completely synchronized state, we present in Fig. 1 the evolution of N = 20 identical Kuramoto oscillators arranged on a wheel graph.A N-vertex wheel graph is formed by N − 1 nodes arranged on a ring and connected to a common vertex.Fig. 1(a) shows the wheel graph representing the connectivity network of the system.The result in Fig. 1(b) displays the evolution of phases θ i (t) for each oscillator i, (i = 1, 2, . . ., 20) as a function of time t, the values of the phases are encoded in the color bar.In this case, the system starts with random initial phases distributed uniformly on the interval [0, 2π) and eventually all of them reach the same value.In Fig. 1(c) we present the behavior of the order parameter r(t), the different curves represent 1000 realizations of the process and the dashed line describes the ensemble average r(t) .Initially, r(t) takes different values between 0 and 0.6 due to the random distribution of the phases at t = 0, then r(t) evolves for all the conditions to 1 as the oscillators reach a coherent state.This general behavior is also observed in the ensemble average r(t) .
Despite the simplicity of the modeling of synchronization with identical oscillators, Eq. (2) includes the topology of the network, so this framework is appropriated to deal with the problem of the influence of topology in the dynamics of coupled oscillators [62,64].Additionally, we consider that the presence of stochastic terms or different frequencies of heterogeneous oscillators can potentially mask the effects of the phenomenon we want to study.In the following, we are interested in the global changes in the synchronization process due to the effects of damage in the links of the coupling network; therefore, this model of identical Kuramoto oscillators results in a convenient starting point.

Synchronization times
The evolution to a completely synchronized state in systems of identical Kuramoto oscillators means that all the oscillators acquire identical phases.Strictly, the time at which these systems achieve a fully synchronized state is infinite; nevertheless, it is possible to relax that condition [67].For the sake of convenience, we analyze the time at which the systems are almost synchronized and the order parameter r(t) takes a determined fixed value r close to 1, we denote this synchronization time as τ 0 .From Eq.
(2), we observe that τ 0 is a variable that depends on the topology of the network and the initial phases of the system, if we consider random phases at t = 0, τ 0 can be seen as an stochastic variable and the statistical analysis of this quantity gives us relevant information of the system.
In the following, we define τ 0 as the time at which the system has reached a state such that the threshold r = 0.99 is crossed from below.To find these times, we perform numerical integration of Eq. ( 2) with random initial phases sampled from a uniform distribution in the interval [0, 2π).We use the fourth order Runge-Kutta method [68] implemented in C++ in the GNU Scientific Library [69].Here, it is important to highlight that several works have shown that in the numerical solution of Eq. ( 2), the results do not depend significantly on the integration method implemented or on the time step ∆t used to obtain numerical solutions [23,70,71].
In Fig. 2, we present the probability density of synchronization times τ 0 for different networks with N = 100 nodes.The probability densities were generated from 10 5 realizations of τ 0 obtained using random initial phases.In Fig. 2(a), we analyze the synchronization times for a particular lattice called knight graph formed by all possible movements of a knight chess piece in a chessboard with dimensions 10 × 10 [72].In Fig. 2(b), we show the results for a wheel graph.Panels 2(c)-(d) show the analysis of two random networks with completely different topologies.Figure 2(c) depicts our findings for an Erdős-Rényi network constructed starting from a set of N = 100 of nodes which are joined by edges whose ends are selected at random with probability p = 0.04 among all the vertices [73], in this model the degree distribution can be approximated by the Poisson distribution [1].Finally, Fig. 2(d) presents the results for a Barabási-Albert network generated using a preferential attachment algorithm so each new node is connected to m = 2 nodes chosen accordingly to their number of connections or degree [74], this algorithm provides an example of the emergence of networks with heavy-tailed degree distributions in terms of the elementary process governing the wiring of new vertices joining the network [1].All the networks explored are plotted as insets and vertical dashed lines show the ensemble average for the values τ 0 .
The probability densities ρ(τ 0 ) show that the τ 0 are close to certain values but the curves are asymmetric.The different curves evidence the influence of the topology of the network in the synchronization process.For each network, we notice that characteristics such as the width and height of the curves are different.It is known that there is an effect  2) considering random initial phases drawn from a uniform distribution in the interval [0, 2π) and evaluating through numerical integration the time τ 0 when the threshold r = 0.99 is crossed from below.Vertical dashed lines represent the ensemble average τ 0 and each network analyzed is presented as an inset.We use bins with size ∆τ 0 = 0.1 in the calculation of the probability densities.
of the structure of the network to achieve synchronization [75,76,77] and the results show that this effect can be observed in the probability densities of τ 0 .Concerning the mean value τ 0 for each network, we notice that this quantity is different for all of them, Fig. 2(a) and Fig. 2(b) exhibit results that differ slightly; however, in Figs.2(c) or (d) the differences are more obvious as a consequence of the relation of τ 0 with the connectivity of the networks [67].These results show that τ 0 is a good candidate to study the performance of synchronization of a network of identical Kuramoto oscillators.In the following part, we will see that this quantity becomes relevant when we explore the synchronization of networks with damage.

Modeling the damage on complex systems
In this section, we present a brief overview of the method developed in Refs.[43,45] to model the process of generation of cumulative damage and aging in networks.In general, the process of accumulation of damage reduces the functionality of the links of the network with fails generated stochastically with a preferential attachment mechanism.The characteristic time scale of this process of accumulation of damage is represented by time T .The time interval at which system receives damage ∆T is long compared to the characteristic time of other dynamical processes that can take place in the system such as a transportation or synchronization process, for ∆T = 1 (measured in units of the system damage accumulation), T = 0, 1, 2, . ... We see for example this difference in the time scales in the activity of a metro system, where the daily activity of users in trains is measured in hours or days but the damage of the infrastructure is only observed at the scale of months or years.We refer the reader to Ref. [78] for a formal mathematical treatment of systems that exhibit multiple time scale dynamics.
The damage in the edge (i, j) of the network is generated at time T with probability π ij (T ) [43] where h ij (T ) is a stochastic integer variable such as h ij (T ) − 1 counts the number of random faults that exist in link (i, j) at time T .Initially, at T = 0, h ij (0) = 1, i.e. there are no faults in any edge in the network.At this point, the algorithm uniformly selects an edge.After damage occurs in any edge, T is increased by 1 and the probabilities of the edges to receive a new fault change according to Eq. ( 4).The entire process generates structures with preferential distribution of damage, where new faults are more likely to accumulate in the most affected edges.Since the damage in edge (i, j) is independent of the damage that receives (j, i), the method generates an asymmetric N × N matrix of weights Ω(T ) describing the couplings between oscillators.Ω(T ) is introduced to characterize the state of damage of the connections of the system and its elements are defined by [43,45] where α ≥ 0 is a real value called the "misrepair" parameter that quantifies the response of the system to damage in terms of repair capacity [43]; in this manner, α shows how detrimental the accumulation of damage is.On the one hand, in the limit α → 0, the system is able to repair completely after receiving damage, in such a way that Ω ij (T ) → A ij .In the limit α → ∞, the system can not be repaired so a hit in a link is equivalent to its removal from the network, this removal of edges gradually brings the system to the percolation limit, where the network becomes disconnected.
This aging process has been used to analyze the transport capacity of networks under stochastic damage and has allowed researchers to conclude that more complex structures are more resilient to the process of aging [43,45].A recent application is the use of this model to study the robustness of metro systems under damage accumulation [46].In the following part, we study the effects of damage accumulation on networks whose main function is to synchronize.In particular, we are interested in cases where the system is able to repair but not in a perfect way so it starts to accumulate damage in its links, which means that α only takes finite values.

Functionality
In Sec.2.1, we introduced the Kuramoto model of identical oscillators in a network without damage.Let us now include the damage in the network and introduce a measure to quantify its consequences in synchronization.We consider that damage modifies the connectivity of the underlying structure, generating a reduction of coupling between the oscillators in the network.In this new structure, the Kuramoto system receives damage in a time interval ∆T which is considerably large in comparison to the time scale t of the synchronization process of the network, ∆T is enough time to the system to synchronize from an incoherent state.The variable T represents the time at the scale of damage generation.But due to the fact that T increases by 1 once the system receives new damage, T also provides information on the amount of global damage that the system has received.The matrix of weights Ω(T ), with elements defined in Eq. ( 5), plays the role of a weighted connectivity matrix in the initial Kuramoto model.In this manner, for each configuration of damage characterized by T , the system of nonlinear equations that define the dynamics of the Kuramoto oscillators is given by Linear models are commonly used to study synchronization.A linear approximation of the Kuramoto model is valid if the system operates close to synchronization.In such scenarios, the evolution of the system can be effectively studied in terms of the eigenvalues and eigenvectors of the Laplacian matrix of the network.To complement this section, the linear approximation is discussed in detail in Appendix 6.1.However, we are interested in the general dynamics described by Eq. ( 6).
Several works have studied the influence of time-varying couplings in the synchronization process of networks in scenarios where the changes of the coupling network can be fast or slow in comparison to the synchronization times of the system [56,57,58,59].In that framework, the effects of the network alteration are evaluated by means of stability analysis of the synchronization state.We implement a different approach assuming the dynamics as a multiple time scale process for which the changes produced by the damage occur at a larger scale in comparison with the synchronization times.
In the following, we assume that the function of the system at each configuration of damage described by Ω(T ) is to reach the global synchronization starting from random phases.In order to evaluate this capacity of synchronization on networks with damage at each T , we propose a measure of "functionality" F s (T ) that quantifies the "health" of the system comparing the synchronization times in Sec.2.2, at different stages of damage.We define Here, we denote as τ (T ) to the synchronization time when the network has suffered damage with couplings described by Ω(T ), the initial phases at t = 0 are chosen randomly but remain the same when evaluating τ 0 and τ (T ).The functionality F s (T ) in Eq. ( 7) is a global measure of the effects of damage in the synchronization process on networks.Particularly F s (0) = 1 and on average decreases with T .In situations when the network fails to synchronize after a certain amount of damage has been added τ (T ) → ∞.Consequently, the network is incapable of carrying out its function leading to F s (T ) → 0.

Effects of damage in synchronization process
In this section, we study the effects of damage on the global synchronization of Kuramoto systems.To illustrate this, in Fig. 3 we present the behavior of synchronization of a wheel graph with N = 20 nodes for different values of T characterizing the damage and α = 0.5.In Fig. 3(a) we depict the network with its weighted links at T = 10 3 .The network is drawn showing only the heaviest edge between each pair of connected nodes, and the weights of the links are codified according to the color bar.We can observe how some of the links keep their original weight while others are notably affected as a consequence of the preferential damage distribution.In Fig. 3(b) we present the evolution of the phases of each oscillator at T = 10 3 .The system starts with the same initial conditions that we used in the example in Fig. 1(b), after some steps the phases reach a coherent state; however, we can observe that in this case, the time to achieve the coherent state is longer than the time used in the network without damage presented in Fig. 1(b).Fig. 3(c) shows the network at T = 10 4 , in this case, we observe a very deteriorated network with connections that almost disappear.Fig. 3 (d) depicts the evolution of phases of the system at T = 10 4 , in a similar way to the case with T = 10 3 , the system starts with the same initial conditions and evolves to a synchronized state, but now the time needed to achieve this configuration is still longer than in the cases with less damage.
On the other hand, Fig. 3(e) displays the ensemble average of the order parameter r(t) for T = 0, T = 10 3 and T = 10 4 .In each case, the average is over 1000 realizations of Monte Carlo simulations, and the set of random initial conditions used to evaluate r(t) in each realization at T = 0 is the same for the cases T = 10 3 and T = 10 4 .The curves show the general behavior of the system under different random initial conditions and diverse stages of damage.The results show that on average, the system evolves from an incoherent initial state to a synchronized state.Initially, all the curves remain close to each other and r(t) does not change significantly, eventually, the curves separate and evolve at different rates so each one reaches 1 at a different time.The curve for T = 0 arrives first, then the curve for T = 10 3 and finally the curve for T = 10 4 .This example shows a system with complete synchronization and how the damage in the coupling network modifies the time in which the system reaches that state of full coherence.

Functionality reduction in the Kuramoto model
In Sec.3.2, we introduced a measure of functionality of the synchronization F s (T ) which gives us insight into the response of the synchronization process under the damage of the network and it can be used to compare networks with different structures.In Fig. 4 we present the evolution of the ensemble average of the functionality F s (T ) defined in Eq. ( 7) as a function of T for the networks explored in Fig. 2 and different values of α.For the evaluation of F s (T ) we use the numerical values of τ 0 and τ (T ) obtained when the system reaches a coherent state with a fixed order parameter r = 0.99 (see the appendix in Sec.6.2 for a detailed discussion of the effect of having different threshold values r).The averages were calculated over 1000 realizations.The different values of α are codified in the color bar that appears in Fig. 4(a).In Fig. 4(a) we present the The Monte Carlo simulations were performed for 1000 realization using random initial phases sampled from a uniform distribution in the interval [0, 2π).
results for the knight graph, in Fig. 4(b) the wheel graph, in Fig. 4(c) the Erdős-Rényi network and in Fig. 4(d) the results for the Barabási-Albert network.In all the cases we observe that F s (T ) is closer to 1 when the damage in the system is small nevertheless, as the damage increases, the functionality decreases.These results are in agreement with previous findings in Fig. 3, showing the increase of the synchronization time τ (T ).Also, the changes in α affect the functionality since the increase in α implies an increase in the intensity of damage so for large values of α the loss of functionality is more significant.One of the interesting features that we find in this result is the similarity observed in Figs.4(b)-(d).Despite the topological differences of the networks, including the different degree distributions, the results of the average ensemble are identical, and the case depicted in Fig. 4(a) is just a re-scaled version.The reason behind this result is the number of edges of the analyzed network.As the damage is distributed across the edges, part of the robustness exhibited by these structures against damage comes from the number of edges.The knight graph has 576 edges, the wheel has 396, and the Erdős-Rényi and Barabási-Albert networks have 392 edges.The ensemble average of the functionality of synchronization scales with the number of edges of the graphs.
In the next section, we continue to study the properties of this measure, especially its relationship with the topology of the networks.

Synchronization on small graphs with damage
Having defined synchronization on networks under the influence of accumulation of damage, as well as introduced the concept of functionality F s (T ) in Eq. ( 7), and its reduction with T ; in this section, we explore the relation between F s (T ) and the structure of the network.To this end, we study the effect of damage on all the connected non-isomorphic graphs with size N = 6, the set of graphs analyzed is available in Ref.
[79] and contains 112 graphs providing a great variety of topologies that includes several trees (for example, the linear graph or the star graph), networks with cycles with different lengths and structures with a high density of edges including the fully connected graph.Also, having only six nodes, we can generate more realizations of the Monte Carlo simulations to perform a better statistical analysis of the results.
In the following, we are interested only in structures that fully synchronize in the absence of damage, thus we omit three particular graphs (a ring with six nodes and two networks formed by a ring with five nodes added one and two edges).In this manner, we consider 109 graphs in the dataset.For each network, we calculate the synchronization times τ 0 and τ (T ), defined as the minimum time t for which r(t) in Eq. ( 3) reach the value r = 0.99 in the synchronization on the original network (to obtain τ 0 ) and the structure with accumulation of damage with T = 100 and α = 0.5 (to obtain τ (T )).We generate 10 7 pairs (τ 0 , τ (T )) using initial random phases.However, we maintain the same initial condition to calculate τ 0 and τ (T ) in each realization.The results are obtained using numerical integration of Eqs.(6).
Once we generated the set of values (τ 0 , τ (T )) for each graph, we can evaluate the ensemble average of the functionality F s (T ) = τ 0 /τ (T ) .Our findings are presented in Fig. 5 where, in panel 5(a), the graphs are sorted in increasing values of F s (T ) ranging from the most affected by damage with F s (T ) = 0.34175 (associated to the star graph) to the structure that better tolerates the damage with F s (T ) = 0.59621 (corresponding to the fully connected graph).The results show variations in the values of the functionality in the set of graphs analyzed, in some cases with differences F s (T ) between two graphs at the order of 10 −5 ; however, such variations can be identified with 10 7 realizations.In addition, in panel 5(b) we plot F s (T ) as a function of the total number of edges |E| ≡ N i=1 N j=1 A ij (considering both directions in each edge).In this visualization of the results, we see that the total number of edges is an important quantity to determine if the network is resilient to cumulative damage; in particular, the higher |E|, the damage is distributed in more edges making them be less affected in its functionality.However, the relation between F s (T ) and |E| is non-linear and, for a fixed value |E|, subtle differences appear associated to the particularities in the topology In Fig. 6, we present all the 109 graphs sorted according to the values of F s (T ) at T = 100, as in Fig. 5.This classification starts with tree networks (1 to 6) with the lowest functionality values showing reduced tolerance to damage in the synchronization process.The following networks are more connected and have structures such as cycles and more redundant paths.Something important to highlight in the results is that this order obtained for synchronization differs from a classification obtained in Ref. [43] for networks with damage but with a global functionality associated with the capacity of random walkers to reach any node of the network.For example, in graph 50 called barbell (formed by two triangles and an edge that joins them), in the transport process, this structure is more fragile since the removal or damage of the edge that joins the two triangles alters drastically the connectivity of two important parts of the network affecting the global transport.In contrast, the implementation of the functionality in Eq. ( 7) considers the couplings between neighbors and global synchronization.Therefore, in comparison with the transport, the synchronization in the barbell graph is not severely affected by the fragility of a link connecting the two triangles, since each group of nodes is independently resistant to damage.This shows that in the effects of accumulated damage the mechanisms that make a network more resilient in transport issues are not the same when the function of the network is to reach states with complete synchronization.
On the other hand, one aspect that stands out among the structures with the same number of edges is the fact that the most resistant topologies have a subgraph with higher connectivity and linear chains of one or two nodes with one of their ends fixed Figure 6: Non-isomorphic connected graphs with N = 6 that synchronize at finite times.The networks are sorted according to their response to damage considering the value of F s (T ) at T = 100 presented in Fig. 5 for α = 0.5 and 10 7 realizations of synchronization on damaged networks.The graphs are obtained from [79].Of the 112 graphs in the original dataset, 3 graphs were discarded because they do not reach complete synchronization for several initial conditions.to the most connected structure.This case is observed, for example, in networks 19 and 34.One possible explanation is that such a mechanism receives little damage in the parts with fewer connections and it is more likely to concentrate damage in the subgraph with higher connectivity.This group of nodes operates as a regulator that prevents from damage key edges whose affectation would change the synchronization times of the whole structure.Finally, in the classification of the networks in Fig. 6, we do not see a particular pattern in the number of cycles with three, four, or more nodes, something that in the case of diffusive transport makes the network more resistant to damage [46].
Our findings for the functionality F s (T ) show that this quantity gives a global characterization of the effect of damage.However, as it is common in the analysis of complex systems, the ensemble average is just a part of the information and it is necessary to introduce a new quantity to see the effect of damage beyond the dependence on the number of edges.To this end, we explore a measure that considers the probability density ρ(τ 0 ) of synchronization times τ 0 and the probability density ρ(τ ⋆ ), where Then, τ ⋆ is a scaled version of τ (T ).In this way, the functionality acts as a factor to put τ 0 and τ ⋆ in the same scale removing the effect of the number of edges |E| in the accumulation of damage.
In the following, we compare the probability densities ρ(τ 0 ) and ρ(τ ⋆ ).To this end we use the Kullback-Leibler divergence, a standard method to calculate the difference between two probability distributions P (z) and Q(z) describing a stochastic variable z [80,81].For continuous distributions, this divergence is given by [80] Here Q acts as a reference distribution.Also, it is important to emphasize that D KL (P ||Q) is not a distance in the sense of a metric since the distance between P and Q is not necessarily the same as between Q and P .Also, from the definition in Eq. ( 9), it is clear that D KL (P ||Q) > 0 and is null when P = Q.The Kullback-Leibler divergence is widely used in data science [82] but also in physical sciences and engineering, for example in the context of turbulence [83], the characterization of dynamical systems [84], or to compare the activity of vehicles in a transportation system [85], just to mention a few examples.
For the statistical analysis of synchronization times we define In Eq. ( 10), (ρ(τ 0 ) + ρ(τ ⋆ ))/2 is the total probability density obtained for the two sets of synchronization times τ 0 and τ ⋆ .In this manner, the reference distribution (ρ(τ 0 ) + ρ(τ ⋆ ))/2 contains total information of the counts of the synchronization times, with this choice we avoid the division by zero in the integral in Eq. ( 9) used to define R KL {ρ(τ 0 ), ρ(τ ⋆ )}.∆τ is the bin size in the calculation of ρ(τ 0 ) and ρ(τ ⋆ ).In this manner, R KL quantifies the variations between ρ(τ 0 ) and the new information associated with the effect of damage.In cases where ρ(τ 0 ) and ρ(τ ⋆ ) are the same, R KL = 0, evidencing that the only effect of damage is the rescaling of synchronization times quantified by F s (T ) .On the other hand, R KL > 0 when ρ(τ 0 ) and ρ(τ ⋆ ) differ, the value R KL increases with the differences between the two probability densities.In Fig. 7 we present the analysis of R KL for the 109 graphs in Figs. 5 and 6.We use the information generated with the 10 7 realizations of times τ 0 and τ (T ) for our model with accumulation of damage with α = 0.5 and T = 100.In Fig. 7(a) we depict the results as a dispersion of points in three dimensions with coordinates (|E|, F s (T ) , R KL ).We analyze the distribution of 109 points using the K-means algorithm to identify a partition of the values into non-overlapping subgroups (clusters).This unsupervised classification method assigns to data points a cluster such that the sum of the squared distance between the data points and the cluster's centroid (arithmetic mean of all the data points that belong to that cluster) is at the minimum [86], the number of clusters K has to be predetermined.We use the Python library Scikit-learn [87] for an implementation of the K-means algorithm.In Fig. 7(a) we show the results obtained for K = 3, we refer to each set of points as groups 1 to 3. The points in each group are presented with different colors.Here it is worth mentioning that other values of K can be explored; however, we center our discussion using three groups.This classification for the data shows graphs in group 1 are affected significantly by damage, having low functionality, and with major changes in the probability densities ρ(τ 0 ) and ρ(τ ⋆ ).In group 2 the F s (T ) increase in comparison to the values in group 1 and the R KL are reduced, revealing that the probabilities densities of synchronization times present moderate variations.Finally, in group 3, the graphs are more resistant to damage, and ρ(τ 0 ) and ρ(τ ⋆ ) suffer minimal variations.To illustrate the features found in the classification, in Figs.7(b)-(g), we discuss the statistical analysis of the synchronization times in three particular graphs, each in one of the groups 1 to 3.
In Figs.7(b)-(c), the linear graph (graph 6 in Fig. 6) is analyzed.The synchronization in this structure is classified in group 1.In panel (b) we show the two-dimensional histogram with the counts of 10 7 pairs (τ 0 , τ (T )), the bin counts are codified in the color bar and the respective network is presented as an inset.The results in this representation show the variations of τ (T ) due to the damage and in general, the relation with τ 0 is nonlinear.In this case, F s (T ) only contains partial information of the relation between τ 0 and τ (T ).Additional modifications associated with the accumulation of damage in the synchronization are observed in panel (c) with the statistical analysis of τ 0 and the scaled time τ ⋆ .Here, it is worth noticing that ρ(τ 0 ) is bimodal, revealing that in the linear graph, some initial conditions generate in high proportion the synchronization with times around τ 0 = 2.5 whereas, from other initial conditions, the system synchronizes around τ 0 = 9.The effect of damage reduces the two relative maximums observed in ρ(τ 0 ).The probability densities differ significantly producing higher values of R KL in Eq. (10).Similar features are found in all the graphs in group 1.
Similarly, in Figs.7(d)-(e), we present the analysis of the graph 40 in Fig. 6 classified in group 2. In panel (d), we see that the pairs (τ 0 , τ (T )) disperse around a line.The analysis of the probability densities ρ(τ 0 ) and ρ(τ ⋆ ) in (e) show moderate differences between the two curves evidenced with an intermediate value of R KL , these characteristics are present in all the graphs in group 2. Finally, we have group 3 in the classification.The effect of damage in this group is illustrated with the analysis of graph 93 in Fig. 6 presented in Figs.7(f) and (g).The results show that the effect of damage is well described by a linear relation, a fact also evidenced in the minimal differences between ρ(τ 0 ) and ρ(τ ⋆ ) producing a R KL closer to zero.
The results in Fig. 7 show that R KL captures the effect of the accumulation of damage in synchronization processes beyond the number of edges.As a final test for this quantity, we evaluate the functionality F s (T ) and R KL for the four networks with N = 100 nodes analyzed in Fig. 2. Our findings are presented in Table 1, where we analyze 10 5 pairs of values generated for T = 1000, α = 0.5 considering synchronization times for r = 0.99.The values of R KL are calculated using Eq. ( 10) with probability densities obtained with bin counts with ∆τ = 0.1, the probability densities for ρ(τ 0 ) were analyzed in Fig. 2.  1: Characterization of synchronization and accumulation of damage for the networks with N = 100 nodes in Fig. 2. The values are obtained from the statistical analysis of 10 5 pairs of synchronization times (τ 0 , τ (T )) generated for T = 1000, α = 0.5 and r = 0.99.

Network
The values reported in Table 1 allow us to classify the networks studied in Fig. 2 in terms of their resistance to damage generated with the preferential attachment algorithm and the modifications in the synchronization times.The networks are sorted according to the values F s (T ) in decreasing order presenting the networks from the most robust to the most fragile.The results show that the knight network is the most resistant with the highest functionality, a fact attributable to the number of edges.Then, we have the Barabási-Albert network with a lower number of edges.In both most resistant networks, we observe small values of R KL showing that in these cases the damage only rescales synchronization times but the probability densities of the scaled time τ ⋆ are similar to the found in Figs.2(a) and (d) for ρ(τ 0 ).This behavior is analogous to the observed in group 3 for synchronization on graphs with N = 6.In contrast, with the same number of edges as the Barabási-Albert network, the Erdős-Rény network has lower functionality but the effect of damage is more clear with the value of R KL that increases.A similar result is found in the wheel.In fact, the Erdős-Rény network and the wheel suffer intermediate modifications in the probability densities of the synchronization times; the characteristics of the changes are similar to those observed in group 2 in the analysis in Fig. 7.

Conclusions
In this paper, we explored the consequences of the accumulation of damage in the synchronization process of networks of identical Kuramoto oscillators.The algorithm implemented generates random cumulative damage in the links producing a directed weighted network modeling the detrimental of the couplings between oscillators.In this system, cumulative damage generates changes in the synchronization times.We propose global measures to quantify the changes in the synchronization process as the damage increases in the network.The first one that we call functionality of synchronization F s (T ) is defined as the quotient between the synchronization time when the network does not have damage and the synchronization time when the network is affected by the damage, in both cases we consider the same initial conditions of the oscillators as random phases.The values of F s (T ) decrease with damage and capture global features of the reduction of the capacity of the networks to reach synchronization; however, when we analyze the ensemble average of this quantity most of the captured information relates to how the number of links in the graph improves its response to damage.The second measure R KL captures the changes that the probability densities of the synchronization times experiment with damage, its value increases when the changes in the probabilities are more noticeable.We apply all those ideas to study the synchronization process in different topologies that admit fully synchronization for most of the initial conditions as a wheel, a knight graph, an Erdős-Rény network, and a Barabási-Albert network.Also, we analyze the set of non-isomorphic connected graphs with six nodes.
The findings of our research show that globally the cumulative random damage re-scales synchronization times of the networks, this effect is measured by means of F s (T ) that depends on a non-linear way of the number of edges of the networks, its ensemble average F s (T ) is able to quantify some differences in the synchronization process of networks with damage beyond the number of edges and due to their topology, nevertheless the deviations are very small and it is necessary many realizations of the process to be able to capture them.On the other hand, the measure R KL uses more information about the synchronization process and evaluates better the effects of cumulative damage.The combination of both measures F s (T ) and R KL results in a classification of networks that synchronize in three groups according to whether their damage tolerance is low, medium or high focusing in their topology rather than their number of edges as the case of the set of non-isomorphic connected graphs with six nodes.
The methods and results of this work contribute to a better understanding of the synchronization, in general, they represent an approach to understanding the degradation of the functions of complex systems that accumulate damage in time.
As a generalization of the model, it would be interesting to consider the process of accumulation of damage commensurable with the synchronization of the system.This problem can be studied in the framework of time-varying coupling networks presented in [56,57,58,59,60] including the analysis of changes of the stability of the synchronization process in networks with damage.Furthermore, since we are studying a model with identical oscillators, new research can include features such as heterogeneous oscillators or the presence of noise.The approach introduced paves the way for the exploration of other dynamical systems in the presence of damage.

Effect of damage on the linear approximation of the Kuramoto model
In the Kuramoto model, it has been shown that for initial conditions closer to the synchronization is valid a linear approximation; in this particular case, synchronization times scale with the inverse of the second non-zero smallest eigenvalue of the Laplacian matrix that defines the connectivity of the network [67,88,89].Also, under this linear approximation of the Kuramoto model of identical oscillators, the stability of the synchronous state can be assessed using a master stability function approach [76,90] that gives the necessary conditions for the system to synchronize and establishes restrictions on the eigenvalues of the Laplacian matrix of the network that guarantee the negativity of the maximum Lyapunov exponent [76,90,91,92].
In this appendix, we present the modifications that introduce the effect of damage in the linear approximation of the Kuramoto model.For small values of θ j (t) − θ i (t), is valid the linear approximation of Eq. ( 6) where the strength S i (T ) of the node i is defined as S i (T ) ≡ N ℓ=1 Ω iℓ (T ) and δ ij denotes the Kronecker delta.In this manner, considering the form of the Laplacian matrix of a weighted network, we define the elements L ij (T ) of the Laplacian matrix L(T ) as Therefore, the linear approximation in Eq. ( 11) defines the dynamical process The integration of Eq. ( 13) leads to Using Dirac's notation for the eigenvectors, we have a set of right eigenvectors {|Ψ j (T ) } N j=1 that satisfy the eigenvalue equation L(T ) |Ψ j (T ) = µ j (T ) |Ψ j (T ) for j = 1, . . ., N. With this information, we define the matrix Q(T ) with elements Q(T ) ij = i|Ψ j (T ) and the diagonal matrix Λ(t, T ) = diag(e −tµ 1 (T ) , e −tµ 2 (T ) , . . ., e −tµ N (T ) ).These matrices satisfy where Q(T ) −1 is the inverse of Q(T ).Using the matrix Q(T ) −1 , we define the set of left eigenvectors { Ψi (T ) } N i=1 with components Ψi (T )|j = (Q(T ) −1 ) ij .Therefore, the solution for the linear dynamics in Eq. ( 14) takes the form Then, the effect of damage accumulation in the linear approximation defines a process that can be explored analytically using the same methods developed in our study of the Kuramoto dynamics; however, all the temporal evolution is determined by the eigenvalues and eigenvectors of L(T ).Here it is worth to notice that the temporal evolution on heterogeneous weighted networks, i.e., in cases where the strength of nodes S i (T ) is not a constant, the linear dynamics is not equivalent to a process of L(T ) in the complex plane for T = 1, . . ., 10 4 codified in the colorbar.We represent with black dots the eigenvalues of the Laplacian matrix without damage L(0).(b) Evolution of the ensemble average of the order parameter r(t) at T = 0 for the network using the linear evolution and the Kuramoto model, average values are obtained considering 1000 realizations.(c) Results obtained for r(t) for the system with damage at T = 1000.(d) Probability densities ρ(τ (T )) of synchronization times τ (T ), the results are generated using 10 6 realizations of the linear and Kuramoto models to obtain the time τ (T ) when the order parameter r(t) in Eq. (3) reaches the value r = 0.99 for T = 0, the same analysis is presented in (e) for T = 1000, bin sizes ∆τ = 0.1 are used.In panels (b)-(d), dashed lines indicate the results for the linear model obtained through numerical evaluation of Eq. ( 16), continuous lines show the results using the numerical integration of the Kuramoto model in Eq. ( 6), both dynamics used random initial conditions uniformly distributed in the interval [0, 2π). of diffusive transport associated to continuous-time random walks defined by the normalized Laplacian with elements L ij (T ) ≡ δ ij − Ω ij (T )/S i (T ) [see Ref. [45] for a detailed discussion of the discrete-time random walk dynamics on networks with accumulation of damage with transition probabilities between nodes i, j given by w i→j (T ) = Ω ij (T )/S i (T )].
To illustrate characteristic features introduced by the effect of damage in the linear dynamics, in Fig. 8, we explore the algorithm for damage accumulation with α = 0.5 on the wheel graph with N = 20 discussed in Fig. 3.In Fig. 8(a), we depict the eigenvalues {µ ℓ (T )} N ℓ=1 of the Laplacian matrix L(T ) for one realization of the algorithm with T = 1, . . ., 10 4 .The sets of eigenvalues are presented in the complex plane with different colors for each T , the eigenvalues for the structure without damage (T = 0) are shown with black dots.In this representation of the results is evident how the asymmetry of Ω(T ) produces complex values of µ ℓ (T ).
On the other hand, in panels 8(b)-(c), we compare the ensemble average r(t) as a function of t for the Kuramoto model in Eq. ( 6) and the linear approximation obtained from the numerical evaluation of the analytical result in Eq. ( 16).The values are obtained for the case without damage T = 0 (panel (b)) and T = 1000 (panel (c)).Using a similar approach, in Figs.8(d The results in Fig. 8 reveal the marked differences between the linear approximation and the Kuramoto model.On average, synchronization times for the linear case are lower than the values observed in the non-linear dynamics.The discrepancies are due to the use of initial random phases for which the approximation in Eq. ( 11) is not valid.However, in other types of problems with initial conditions closer to synchronized phases, the linear dynamics with phases in Eq. ( 16) gives a good approximation for the synchronization where we can see directly that the eigenvalues of the Laplacian matrix are associated with characteristic times of the process.In our general modeling, we are interested in a functionality measure defined in terms of the capacity of a system to reach synchronization from random initial phases.To this end, in the discussion presented in the main text, we restrict our study to the Kuramoto model.

Variations of the functionality with the threshold value r
In our study of the effects of damage in synchronization in the main text, we selected a particular threshold value of the order parameter r in Eq. (3).It was convenient to choose the threshold r = 0.99 to study configurations close to complete synchronization (obtained when r = 1).Nevertheless, it is important to understand the effects of choosing different values of r and how the threshold selected can affect our findings.In this appendix, we explore numerically the implications of the threshold value r in the results of the ensemble average of the functionality F s (T ) defined in Eq. (7).We explore two networks for two values of the parameter α and 1 ≤ T ≤ 10 4 .
We define ratio η r (T ) ≡ F s (T ) r F s (T ) r=0.99 , where F s (T ) r denotes the ensemble average of functionality for a given threshold value r; in particular F s (T ) r=0.99 is obtained with r = 0.99, the value we used in the different analyses discussed in the main text.In Fig. 9, we present numerical results obtained for η r (T ) using Monte Carlo simulations, the methods implemented are the same as described for Figs. 2 and 3 but now using the threshold r = 0.6, 0.7, 0.8, 0.9 and 0.99 to define the synchronization times and α = 0.5, 1 for the wheel graph and the Barabási-Albert network with size N = 100 studied in Fig. 2, the ensemble averages were obtained using 1000 realizations.
In the panels in Fig. 9, we plot η r (T ) as a function of T , each distribution of points is obtained with a particular value of r codified in the markers in panel (a).The dots depict the case η r (T ) = 1 when F s (T ) r coincides F s (T ) r=0.99 .In Figs.9(a)-(b) we depict the results for the wheel graph using α = 0.5 and α = 1.0, respectively.Similarly in Figs.9(c)-(d), we present the results for the Barabási-Albert network.In both cases, the numerical values reveal that for α = 0.50, η r (T ) are close 1 showing to F s (T ) r are similar for the different values of r to the result obtained for r = 0.99.This is evident in both the structures, the wheel and the Barabási-Albert network.However, in panels (b) and (d) generated for α = 1, the changes of F s (T ) r are more noticeable evidenced by higher values of η r (T ) for T ≥ 10 3 .The results in Fig. 9 characterize the effect of the threshold r and how affects the values of F s (T ) and show that for r close to 1.0, maintaining α < 1, the F s (T ) are not significantly affected by the choice of r.

Figure 1 :
Figure 1: Synchronization processes in a wheel graph with N = 20 nodes.(a) Wheel graph.(b) Temporal evolution of the phases θ i (t) as a result of the numerical integration of Eq. (2) using random initial phases sampled from a uniform distribution in the interval [0, 2π), the value of the phases is codified in the color bar.(c) Evolution of the order parameter r(t) in the synchronization of the network.Thin lines represent the results of 1000 realizations using random initial phases and the dashed line shows the ensemble average denoted as r(t) .

Figure 2 :
Figure 2: Probability densities ρ(τ 0 ) of synchronization times τ 0 for networks with N = 100 nodes.(a) A knight network with dimensions 10 × 10, (b) a wheel, (c) an Erdős-Rényi network with p = 0.04, (d) a Barabási-Albert network with m = 2.The results are generated using 10 5 Monte Carlo simulations of the model in Eq. (2) considering random initial phases drawn from a uniform distribution in the interval [0, 2π) and evaluating through numerical integration the time τ 0 when the threshold r = 0.99 is crossed from below.Vertical dashed lines represent the ensemble average τ 0 and each network analyzed is presented as an inset.We use bins with size ∆τ 0 = 0.1 in the calculation of the probability densities.

Figure 3 :
Figure 3: Synchronization process in a wheel with N = 20 nodes with damage at T = 10 3 and T = 10 4 with α = 0.5.Panel (a) shows the network and its weighted connections at T = 10 3 , only the links with the largest weight are presented, values are codified in the color bar.Panel (b) depicts the evolution of the phases as a function of t at T = 10 3 .The color represents the value of the phases θ i (t) codified in the color bar.Panels (c) and (d) present the state of the network at T = 10 4 and its evolution.Panel (e) displays the evolution of the ensemble average of the order parameter r(t) for the network at the different scenarios of damage, average values are obtained considering 1000 realizations.

Figure 4 :
Figure 4: Reduction of the ensemble average of the functionality of synchronization F s (T ) as the damage increases in networks of different topologies with N = 100.The curves correspond to F s (T ) for 1000 realizations and different values of α encoded in the color bar.We analyze the networks in Fig. 2: (a) Knight graph of dimension 10 × 10, (b) a wheel, (c) an Erdős-Rényi network and, (d) a Barabási-Albert network.The Monte Carlo simulations were performed for 1000 realization using random initial phases sampled from a uniform distribution in the interval [0, 2π).

Figure 5 :
Figure 5: Ensemble average functionality F s (T ) for cumulative damage with α = 0.5 and T = 100 for 109 non-isomorphic connected graphs with N = 6 nodes.(a) Graphs sorted with the values of F s (T ) .(b) F s (T ) as a function of the number of directed edges |E|.The values are obtained considering 10 7 realizations.See details in the main text.

Figure 7 :
Figure 7: Comparison of synchronization on networks with N = 6 under the influence of damage accumulation with T = 100 and α = 0.5.(a) Dispersion of values (|E|, F s (T ) , R KL ) for the 109 networks presented in Fig. 6.Numerical values of R KL for each graph are obtained using Eq.(10) to compare the probability density of the times ρ(τ ⋆ ) of the scaled synchronization time τ ⋆ in Eq. (8) for the network with damage and the probability density ρ(τ 0 ) for the synchronization time on the original structure.(b) Two-dimensional histogram for 10 7 pairs of values (τ 0 , τ (T )) for a linear graph presented as inset, frequency counts are codified in the color bar.(c) Probability densities ρ(z) for synchronization times z = τ 0 and z = τ ⋆ , we use bins with size ∆τ = 0.1.Panels (d) and (e) show the results for graph number 40 and (f) and (g) for graph 93 in Fig. 6.The groups in (a) are obtained using the K-means algorithm (see the main text for a detailed description).

Figure 8 :
Figure 8: Effect of damage with α = 0.5 for the linear approximation and the Kuramoto dynamics on a wheel with N = 20 nodes.(a) Set of eigenvalues {µ ℓ (T )} Nℓ=1 of L(T ) in the complex plane for T = 1, . . ., 10 4 codified in the colorbar.We represent with black dots the eigenvalues of the Laplacian matrix without damage L(0).(b) Evolution of the ensemble average of the order parameter r(t) at T = 0 for the network using the linear evolution and the Kuramoto model, average values are obtained considering 1000 realizations.(c) Results obtained for r(t) for the system with damage at T = 1000.(d) Probability densities ρ(τ (T )) of synchronization times τ (T ), the results are generated using 10 6 realizations of the linear and Kuramoto models to obtain the time τ (T ) when the order parameter r(t) in Eq. (3) reaches the value r = 0.99 for T = 0, the same analysis is presented in (e) for T = 1000, bin sizes ∆τ = 0.1 are used.In panels (b)-(d), dashed lines indicate the results for the linear model obtained through numerical evaluation of Eq. (16), continuous lines show the results using the numerical integration of the Kuramoto model in Eq. (6), both dynamics used random initial conditions uniformly distributed in the interval [0, 2π).
)-(e), we present the statistical analysis of the time τ (T ) necessary to reach for the first time the order parameter r = 0.99 for T = 0 (panel (d)) and T = 1000 (panel (e)).In all the realizations in Figs.8(b)-(d)we consider random initial conditions for the phases θ i (0), with values chosen random (uniformly distributed) in the interval [0, 2π).

Figure 9 :
Figure 9: η r (T ) as a function of T using the threshold r for a wheel graph and a Barabási-Albert network with N = 100.The ensemble averages in η r (T ) defined by Eq. (17) are generated using 1000 realizations for F s (T ) r with threshold values r = 0.6, 0.7, 0.8, 0.9 and compared with the result obtained for r = 0.99 for 1 ≤ T ≤ 10 4 .Wheel graph with (a) α = 0.5 and (b) α = 1.0.Barabási-Albert network with (c) α = 0.5 and (d) α = 1.0.