Non-monotonic transients to synchrony in Kuramoto networks and electrochemical oscillators

We performed numerical simulations with the Kuramoto model and experiments with oscillatory nickel electrodissolution to explore the dynamical features of the transients from random initial conditions to a fully synchronized (one-cluster) state. The numerical simulations revealed that certain networks (e.g., globally coupled or dense Erdős–Rényi random networks) showed relatively simple behavior with monotonic increase of the Kuramoto order parameter from the random initial condition to the fully synchronized state and that the transient times exhibited a unimodal distribution. However, some modular networks with bridge elements were identified which exhibited non-monotonic variation of the order parameter with local maximum and/or minimum. In these networks, the histogram of the transients times became bimodal and the mean transient time scaled well with inverse of the magnitude of the second largest eigenvalue of the network Laplacian matrix. The non-monotonic transients increase the relative standard deviations from about 0.3 to 0.5, i.e., the transient times became more diverse. The non-monotonic transients are related to generation of phase patterns where the modules are synchronized but approximately anti-phase to each other. The predictions of the numerical simulations were demonstrated in a population of coupled oscillatory electrochemical reactions in global, modular, and irregular tree networks. The findings clarify the role of network structure in generation of complex transients that can, for example, play a role in intermittent desynchronization of the circadian clock due to external cues or in deep brain stimulations where long transients are required after a desynchronization stimulus.

Here the state of each oscillator is represented by a scalar variable θ n , its phase, while the interaction between oscillators is described by a square adjacency matrix A nm and a coupling strength κ. We assume that A nm = 1 if there is a link between the nth and mth oscillators, and A nm = 0 otherwise. All the links are supposed to be undirected, therefore matrix A nm is symmetric. Then for each oscillator we define its degree by the formula e iθ n (t) .
A more detailed information is obtained in the following way. We fix the network and run 10 4    from a small region close to zero in the panel (f)) look as monotonic functions (or nearly monotonic functions with only small local maxima/minima); therefore we call these transients to synchrony monotonic.
It turns out that monotonic transients to synchrony are not generic for all types of networks. If we consider a modular network consisting of two equal size all-to-all coupled populations connected by a bridge of two oscillators, figure 2(a), then we find that the system approaches synchrony in a qualitatively different way. In this case the histogram of transient times is bimodal (see figure 2(b)). It has a smaller and narrower peak at τ ≈ 1 and a larger and broader peak at τ ≈ 33. Moreover, if we consider the longest trajectory of the order parameter R(t) we find that it is not a monotonic function. More precisely, we observe the following scenario. The order parameter first grows quickly to a relatively large value R ≈ 0.75, then it decreases back to zero and only afterward it turns to increase again and approaches eventually the asymptotic value R = 1. (For comparison, as shown in figure 2(c), the typical trajectory with mean transient time is monotonic.) The snapshots in figure 2(d)-(g) explain the microscopic mechanism of this non-monotonic transient to synchrony. At the first maximum of R(t) both populations are almost synchronized (recall that the phases −π and π are identical), but the two oscillators on the bridge remain distant from the synchronous cluster (d).
Instead of attracting these two oscillators the phases of two populations start to drift apart from each other (e) and approach an anti-phase configuration (f). Only after such a detour all oscillator phases converge to the completely synchronous state (g). Motivated by the above example we adopted the following definition of non-monotonic transients. Given a system trajectory with R(t) → 1 for t → ∞, one has to find its largest local maximum R max = R(t max ) and its smallest local minimum R min = R(t min ). If the both extrema exist and R max − R min ΔR = 0.2 and t min > t max , then the trajectory is called non-monotonic transient to synchrony. Otherwise it is called monotonic transient to synchrony. Using this definition we classified all transients in the histogram in figure 2(b). We found that the probability of observing non-monotonic transients correlates with the transient length τ . In particular, there exists a critical value τ c ≈ 33.7 such that all transients with τ < τ c are monotonic. On the other hand, non-monotonic transients can be found for τ > τ c only and the chances of observing them increase with τ . All transients with τ close enough to the maximal transient length turned out to be nonmonotonic. Comparing different non-monotonic transients we also noticed that the longer transient is, the more pronounced is its oscillation, see figure 2(c).

Robustness of non-monotonic transients to synchrony
Looking at the network in figure 2(a) one can think that the non-monotonic transients are relevant to its specific symmetries. In order to clarify this issue we repeated our numerical simulations for several modified networks. We kept the system size N = 20 constant and varied the network architecture. First, we doubled the number of oscillators on the bridge connecting two populations (respectively, we decreased the size of each population by one). We found that the non-monotonic transients to synchrony are still observed, see figures 3(a)-(c). Moreover, the distribution of their transient times is described by a similar bimodal histogram but with a larger mean transient time. Next, we generated a random network removing half of the links in each population of the network shown in figure 2(a). We found that the non-monotonic transients persist in the new system as well, see figures 3(d)-(f). Finally, we considered a modular network with populations of different size, figure 3(g), and found again that this modification has no significant impact on the qualitative features of non-monotonic transients to synchrony, figures 3(h) and (i).
Summarizing these results we conjecture that non-monotonic transients to synchrony are typical for networks consisting of two oscillator populations densely connected inside each of them and mutually connected by a chain of other oscillators. On the other hand, it seems that the phenomenon of non-monotonic transients is suppressed in the case of two or more interconnecting chains. For example, figure 3(k) shows a unimodal histogram of transient times obtained for a multiplex network. Note that in this case, only a few of 10 4 random initial conditions result in non-monotonic transients that are exemplified in panel (l).
Apart from the network architecture the occurrence of non-monotonic transients can also be influenced by the oscillators heterogeneity and the network size. Because of the deterministic nature of system (1), small disturbances of intrinsic frequencies ω n cause only minor changes of the order parameter graphs and hence of the transient time histograms. On the other hand, if the frequency variation is of the same order or larger than the coupling strength κ, then the network dynamics can change significantly. In particular, the system can become multistable such that for different initial conditions it converges to different phase locked states. Moreover, the long term dynamics of the order parameter R(t) can become periodic, quasiperiodic or even chaotic what blurs the difference between non-monotonic and monotonic transients. All these issues go beyond the scope of the present paper and we leave them to future work.
Regarding the network size, we tested its impact in the following way. We took the random network in figure 3(d) and constructed analogs for N = 40, 60, 80 and 100, increasing the population size and the number of bridge elements proportionally to N. For each of the new networks we computed its transient time histogram distinguishing between monotonic and non-monotonic transients according to the definition from section 2.1. The obtained results are shown in figure 4. They illustrate two qualitative trends. First, the relative frequency of non-monotonic transients increases for increasing N. Second, the non-monotonic behavior of R(t) becomes more sophisticated for larger N. In particular, it can include more than one oscillations between the instances of higher and lower synchrony, see figure 4(h).

Laplacian matrix eigenvalues and eigenvectors
Considering the stability of the completely synchronous state one often refers to the spectrum of the Laplacian matrix [27,29,37], which is defined by the expression 3 L nm = A nm − k n δ nm , 3 Note that in some textbooks the definition of the Laplacian matrix has the opposite sign. where k n is the degree of the nth oscillator and δ nm is the Kronecker delta. The Laplacian matrix is negative semidefinite, therefore all its eigenvalues λ n are real and nonpositive. Usually one orders them in the decreasing order Then for a connected network, the largest eigenvalue λ 1 vanishes, while all other eigenvalues are strictly negative. Linearizing system (1) in the vicinity of the completely synchronous state one obtains a linear system for small phase perturbations φ n . This system implies that after the aggregation of the phases, the relaxation to synchrony occurs along the slowest stable direction where (φ (2) 1 , . . . , φ (2) N ) T is the eigenvector of L nm corresponding to the eigenvalue λ 2 . For the considered networks where non-monotonic transients to synchrony were observed, there was a pronounced spectral gap between the second and the third eigenvalues, λ 3 /λ 2 1, of the Laplacian matrices (see figure 5). In addition, the eigenvector elements (φ (2) 1 , . . . , φ (2) N ) T resembled the phase snapshots in the model (1) en route to the complete synchronization, compare figures 2(d)-(g) and the insert panel in figure 5(a). These properties of the Laplacian matrix (spectral gap and modular structure of the second eigenvector) seem necessary for the observation of non-monotonic transients in model (1).
Remarkably, the characteristic time scale 1/|λ 2 | is relevant not only to the local dynamics of model (1) in the vicinity of the completely synchronous state, but also to the global dynamics of this model. In figure 6(a) we plotted |λ 2 | τ for all numerical simulations reported above, and found that the networks with the  Note that networks without spectral gap (e.g., globally coupled network), do not follow the relationship. We can understand this different behavior using the following heuristic argument. The transient time is the sum of the aggregation time by which the phases come together, and a relaxation time governed by the dominant eigenvector and eigenvalue of the system. For networks with spectral gap, the transient times are dominated by the relaxation time, and thus the scaling is recovered. In the other networks (e.g., global coupling) the aggregation and relaxation process cannot be separated so easily, and thus deviations occur. Nonetheless, the aggregation and relaxation processes will need future investigations to better understand the dependence of the mean transient times on the network structure. We also computed the standard deviations σ of transient times for all histograms in figures 1-4, 13 and found that the ratios σ/ τ lie in the interval [0.3, 0.5], figure 6(b). In particular, for all cases of non-monotonic transients to synchrony we obtained σ/ τ ≈ 0.5.

Impact of initial conditions on the non-monotonic transients to synchrony
In the experiments described later, because of technical restrictions, we could not generate initial conditions corresponding to randomly distributed phases in this system. However, it was relatively easy to obtain initial conditions, which are random permutations of a completely desynchronized state where the elements populate the limit cycle uniformly with a constant phase difference. Motivated by this fact we checked if the statistics of transient times in the Kuramoto networks is sensitive to such a limitation on the initial conditions. For this, we repeated 10 4 numerical simulations for the model (1) with the modular network from figure 2(a). In the simulations we used the following initial conditions: each time we generated a random permutation of N = 20 elements using the Fisher-Yates shuffle algorithm and reordered the indices of the oscillators in the desynchronized state according to the obtained permutation, figure 7(a). It turned out that the transient time histogram for the new initial conditions, figure 7(b), differs qualitatively from the histogram in figure 2(b). In particular, it is no more bimodal but has a single left-skewed peak for τ ≈ 36. On the other hand, the order parameter graph of the longest detected transient trajectory for the new initial conditions, figure 7(c), is very similar to that one shown in figure 2(c). This allows us to conclude that these initial conditions are sufficient to reveal non-monotonic transients to synchrony in a model (1), but would not reveal the bimodal structure of the corresponding transient time histogram.

Experimental procedures
The experimental setup is a three-electrode electrochemical cell with an electrode array of twelve 1 mm diameter nickel wires as working electrodes, a Hg/Hg 2 SO 4 /sat. K 2 SO 4 reference (RE), and a Pt coated Ti rod as counter (CE) electrode; see figure 8(a) for a schematic of the experiments. The electrolyte is 3 M sulfuric acid solution at 10 • C. The nickel wires were embedded in epoxy, so that only the ends of the wires were exposed to the electrolyte solution. At a constant circuit potential and an external resistance (R ind = 1 kohm) attached to each nickel wire, the electrodes exhibit periodic current oscillations that arise through a Hopf bifurcation [38]. The current of each electrode was measured and digitized from the potential drops across the R ind at a rate of 200 Hz with an NI PXI 6255 data acquisition board. The circuit potential was set to 40 mV above the Hopf point; because of the somewhat different experimental conditions (e.g., slightly different surface areas, changing surface properties of the oxide films), the Hopf point varied in the range V = 1060-1075 mV (with respect to RE) throughout the 200 experiments reported in this work. Under these conditions, for a typical experimental data corresponding to uncoupled oscillators, the mean natural frequency was ω = 0.419 Hz and the standard deviation (for the twelve oscillators) was 3 mHz.
A given network was obtained by adding cross resistors (R c ) between the coupled electrode pairs. The cross resistors induce an electric coupling, so that current can flow between wires due to difference in the electrode potentials; the coupling strength (conductance) is K = 1/R c [26,39]. Three different coupling strengths were used with R c = 1 kohm, 2.5 kohm or 7.5 kohm. Five different networks (each with twelve electrodes) were generated: globally coupled network, modular networks with one bridge interconnection of lengths 3 (figure 8(b)) and 4 (figure 8(c)), a modular network with 4 bridges with one node (figure 8(d)), and an irregular tree network ( figure 8(e)).
Before the coupling was turned on, the oscillators were desynchronized by perturbations of the circuit potential (δV) proportional to the total current (I): where L is the feedback gain and o is an offset calculated as the mean total current. As it was shown in previous studies [40,41], such feedback induces a phase repulsive global coupling in the population and results in a completely desynchronized state, where the phases of the oscillators are nearly evenly distributed. After the desynchronized state was achieved with L = 10 mV mA −1 , the feedback was turned off, and the coupling was turned on after a 5 s wait time. After the coupling was turned on, the currents were recorded until full synchronization was observed. Typical length of the recorded time series is 50-100 s, but full synchronization was achieved in all experiments within 300 s. The phases of the oscillators were calculated with the peak finding method with linear interpolation between the peaks [3].

Global coupling: monotonic transient to synchrony
The synchronization transient with global coupling is shown in figure 9 at a coupling strength K = 133 μS. Figure 9(a) shows the time series of the current with a zoom to the early state of the resynchronization process (first 10 s); the corresponding space-time plot of the current is in figure 9(b). When the coupling was turned on (t = 0 s), the oscillations become quickly synchronized and the currents nearly overlapped by t = 5 s (or two oscillatory cycles). In the early stage of the synchronization, the amplitudes of the oscillations greatly diminished. Between 5 s < t < 25 s, the oscillations remained synchronized, and the amplitudes slowly increased. After about t = 25 s, the system reached full synchrony, where the amplitudes in the synchronized state were similar to those without coupling. As the R vs t plot shows in figure Figure 10 shows  figure 10(e)). At t = 0 s, figure 10(c), the system is desynchronized. In about 2 oscillation cycles, t = 5 s, the order parameter reaches a local maximum, and the elements in modules synchronized in phase, but have anti-phase synchronization between modules; the bridge elements at this stage have phases similar to those in module 2. At t = 12 s (at the minimum of the R vs t plot, figure 10(f)), the two modules remained synchronized in anti-phase configuration, but the bridge element oscillator 5 connected to module 1 now established a phase close to module 1. After this stage, the order parameter increased; at t = 50 s ( figure 10(g)), around halfway to full synchrony, the elements in the modules remain in-phase synchronized, but the phase difference is less than π, and the bridge elements smoothly connect the phases of modules (similarly to a standing wave oscillation). As the full synchrony is achieved at τ = 108.4 s, the phase difference between the oscillations in the modules further decreased, which can be in the space-time plot of the currents in figure 10(a).

Modular network with bridge elements: non-monotonic transient to synchrony
The histogram of the recorded transient times τ in 31 trials is shown in figure 10(h). The mean transient time was found to be τ = 100.4 s with a relative standard deviation of about 0.60. While we still observed a unimodal distribution, there seems to be a long tail to the large transient times.

Comparisons of transients in different networks
To further explore the characteristics of the transients to synchrony, we performed four other sets of experiments with different network topologies with a coupling strength K = 1000 μS. (Note that the coupling strength was stronger in these experiments than those discussed above.) The four network topologies are shown in figures 8(b)-(e). Figure 11 shows the histograms (right column) of the transient times τ , and the R vs t graph for the longest recorded transient (left column). For the two-module network with three bridging nodes ( figure 8(b), same network as discussed in figure 10), the transient time histogram (figure 11(e)) shows that the mean time τ decreased to 16.1 s (from 100.4 s with K = 133 μS). Note that it was expected that the transient time decreased due to the increased coupling strength, and the extent of decrease in the transient time, 6.2 times, is similar to the 7.5 fold increase in the coupling strength. However, at this increased coupling strength, the order parameter did not show clear non-monotonic transient to synchrony in figure 11(a). This observation shows that at this increased coupling strength, the amplitude dynamics could play a role in the transient dynamics.
In search of non-monotonic transient to synchrony at this stronger coupling, the experiments were repeated with a network with four bridge elements; because the total number of oscillators was limited to 12, this network consists of 4 elements in the two modules. A movement of elements from a module to the bridge elements had large impact on the transient characteristics: the mean transient time increased to τ = 22.6 s (see figure 11(f) for the histogram) and the R vs t graph for the longest transient (τ = 41.6 s) exhibited a minimum at t = 7.9 s. We thus see that addition of bridge elements can generate longer transients and give rise to non-monotonic synchronization transients. For further clarification of the role of the bridge elements, four-bridge network was reorganized such that there are four single bridge elements that connect the two modules ( figure 8(d)). In this setup the synchronization transients were very fast ( τ = 4.9 s, see figure 11(g) for the histogram) and the R vs t graph showed monotonic transient to synchrony (figure 11(c)). Similarly, to show the effect of the modular network topology, the experiments were also performed for an irregular tree network (figure 8(e)). We found that τ = 12.9 s and the transients did not generate large desynchronization peaks (see figures 11(d) and (h)).

Network structure and transient time statistics
The six set of experiments presented in figures 9, 10, and 11 were augmented with two additional sets: a global coupling at K = 1000 μS, and the modular network with a bridge of three elements at K = 400 μS. For these eight set of experiments the theoretical predictions on the correlations between the mean transients times and the relative standard deviations with the inverse of the absolute value of the second largest eigenvalue can be experimentally tested. Because the experiments were performed with three different coupling strengths, the mean transient time was rescaled with K/K s , where K is the coupling strength at which the experiment was performed, and K s = 1000 μS. This rescaling reflects that at stronger coupling the transient times are expected to be weaker. Figure 12(a) shows that the product of the rescaled mean transient times τ and |λ 2 | remains nearly constant for the networks with spectral gaps, but there are large deviations with global coupling. Figure 12(b) shows the relative standard deviations; because the standard deviation is relative to the mean transient time, no rescaling was necessary. As predicted by the simulations, with the globally coupled networks the relative standard deviations are smaller (0.13-0.25) than with the other networks where the points fall in the region between 0.4 and 0.6.

Discussions and conclusions
The main dynamical characteristics of the transients to synchronization were studied in Kuramoto networks and electrochemical oscillators.
Numerical simulations with the Kuramoto oscillators identified networks that exhibit relatively simple, predominantly monotonic transients with nearly unimodal transient time histograms. In contrast, certain modular networks connected with bridge elements, can exhibit non-monotonic transients where the Kuramoto order parameter can exhibit a local minimum en route to full synchronization. This minimum implies that during the transient to full synchronization, desynchronization steps that diminish the amplitude of collective oscillations could occur. Microscopically, the mechanism of the desynchronization is the generation of synchronized modules that can take up anti-phase configurations before transition to full synchronization.
In these modular networks, the transient time histograms appeared to be bimodal for random initial conditions. In contrast, the bimodal character of the histogram disappeared with a uniformly distributed initial condition, while the non-monotonic transient prevailed. The mean transient times exhibited a linear relationship with the inverse of the modulus of the second largest eigenvalue of the network. While such relationship would be certainly expected for relatively simple transients, it is of note that complex transients did not fundamentally break this relationship. The relative standard deviations of the transient times were found to be The experiments with electrochemical oscillators confirmed many predictions of the simulations. One experimental challenge was the generation of random initial conditions from a given statistical distribution. The initial conditions were generated with an external feedback of the uncoupled system, that generated a desynchronized initial condition in which the elements had nearly uniformly distributed phases. With global coupling relatively fast transients were seen with monotonic transients. However, modular networks coupled with a single bridge exhibited non-monotonic transients that occurred through intermittent anti-phase synchronization of the modules. In the experiments strong coupling effects could also be observed; for example, at very strong coupling the non-monotonic transients disappeared. The experiments also confirmed that for the network with spectral gaps the mean transient times correlate well with |λ 2 | −1 and that the non-monotonic transients tend to increase the standard deviation of the transient time.
It is noteworthy to compare the results to a previous study [31], where very long transients were obtained in a ring of electrochemical oscillators with a few short cuts similar to a small-world network. In that work, before the (fully) synchronized state was achieved, the system itinerated over various unstable rotating wave states; for each state the order parameter was low representing a nearly desynchronized state in terms of the Kuramoto order parameter. After a few of such phase reordering, the system reached the basin of attraction of the one-cluster state, and the system very quickly converged to the final state. In other words, the aggregation of the phases was a long process (which determined the long transient times), which was followed by a short relaxation. In the transients presented here, there are no rotating waves, the aggregation time is relatively short, and there is a long relaxation time. We thus see that there are probably multiple mechanisms that can generate complex transients in networks of oscillatory processes.
The simulations and experiments presented here can provide a stepping stone for further explorations of the dynamical features of the transients. Some prominent open problems include the exploration of the effect of system size on the features, the theoretical characterization of the transients (e.g., prediction of mean transient times and their standard deviations) and the determination of network properties that are essential for the non-monotonic transients. For example, non-monotonic transients occur not only in two population modular networks, but also in networks composed of three populations, see figure 13. It seems likely that they can also be found in the case of more populations, though their properties can vary from case to case. Another open problem concerns the fact that in many (but not all) of the networks, in which non-monotonic transients were found, there was a spectral gap between second and third eigenvalues of the network (λ 3 /λ 2 1). In such networks, large rearrangements of the initial conditions could be required to reach the final stage of resynchronization, where the phase should be aligned along with the second eigenvector of the network Laplacian matrix. Experimental identification of non-monotonic transients could aid the location of network modules from dynamical measurements where the network structure is not known, for example, through network reconstruction techniques [42].

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.