How clock heterogeneity affects synchronization and can enhance stability

The production process of integrated electronic circuitry inherently leads to large heterogeneities on the component level. For electronic clock networks this implies detuned intrinsic frequencies and differences in coupling strength and the characteristic time-delays associated with signal transmission, processing and feedback. Using a phase-model description, we study the effects of such component heterogeneity on the dynamical properties of synchronization in networks of mutually delay-coupled Kuramoto oscillators. We test the theory against experimental results and circuit-level simulations in a prototype system of mutually delay-coupled electronic clocks, so called phase-locked loops. Contrary to the hindering effects of component heterogeneity for the synchronization in hierarchical networks, we show that clock heterogeneities can enhance self-organized synchronization in networks with flat hierarchy. That means that beyond the optimizations that can be achieved by tuning homogeneous coupling strengths, time-delays and loop-filter cut-off frequencies, heterogeneities in these system parameters enable much better optimization of perturbation decay rates, the stabilization of synchronous states and the tuning of phase-differences between the clocks. Our theory enables the design of custom-fit synchronization layers according to the specific requirements and properties of electronic systems, such as operational frequencies, phase-relations and e.g. transmission-delays. These results are not restricted to electronic systems, as signal transmission, processing and feedback delays are common to networks of spatially distributed and coupled autonomous oscillators.


I. INTRODUCTION
In modern electronics the coordination of complex systems and processes in time is necessary for a defined system behavior, efficient operation and reliable parallel information processing [1][2][3]. To establish such coordination via a global time reference for spatially distributed clock elements can be a challenging task, especially at high frequencies and in the presence of signal transmission-and processing-delays. Such delays are induced by finite signal propagation speeds in transmission lines and due to signal processing, e.g., filtering. Furthermore, integrated electronic systems display a considerable degree of heterogeneity in their component characteristics due to, e.g., the production process of semiconductor technology [4][5][6]. This has far reaching consequences for the architectures of electronic systems, requiring their functionality to be robust against such heterogeneity. Example systems are global and indoor positioning, large antenna, radar and sensory arrays, multi-processor computer architectures, terahertz based technology and databases on the internet [1]. A common approach to synchronization in such systems is to entrain imprecise electronic clocks, so called phaselocked loops (PLLs) [7], hierarchically with a dedicated and precise reference clock (usually a quartz). Such references feed their signal unidirectionally into a clocktree that becomes increasingly complicated as the system size grows [8,9]. Therefore such systems are often synchronized only locally by globally asynchronous, locally synchronous operations (GALS ) [10]. A novel approach to the synchronization of large spatially distributed elec-tronic systems is to allow the formation of self-organized synchronous states [11][12][13][14]. This is inspired by robust self-organized synchronization without hierarchical structures as found in biological systems, where synchronization is achieved robustly in highly noisy environments with strong heterogeneities and in the presence of considerable time-delays [15,16]. So far, such networks have been studied analytically for homogeneous clock networks and and have been tested experimentally with weakly heterogeneous electronic prototype clocks [5,13,17]. To enable this technology for applications in electrical engineering, it is important to understand the consequences of heterogeneity on the collective selforganized dynamics in finite-size systems [14]. Furthermore, studying synchronization of heterogeneous clocks is relevant beyond electronic systems and can provide insight in, e.g. power grid and biological neural networks, consisting of strongly heterogeneous units [19][20][21]. In this paper we study and analyze the effects of heterogeneity on the synchronization dynamics of mutually delay-coupled electronic clocks. We use a Kuramoto-type model, i.e., networks of delay-coupled phase oscillators with node dynamics that include processing and inversion of signals and delayed feedback [22,23].
The paper is organized as follows. We introduce a model for heterogeneous networks of non-identical mutually delay-coupled digital electronic clocks with signal filtering in Sec. II. For such systems we calculate the frequencies and phase configurations of synchronized states, and analyze their stability as a function of different heterogeneous system parameters in Sec. III. In Sec. IV we discuss the effects of these heterogeneities and compare to arXiv:1908.11085v2 [nlin.AO] 2 Sep 2019 the case of networks of identical clocks. We present how the different types of time-delays and their interplay affects the dynamical properties of synchronization in such systems. We summarize in Sec. VI, connect the results to modern electronic applications and components, and discuss the potential of self-organized synchronization for microelectronic systems and processes.

II. MUTUALLY DELAY-COUPLED ELECTRONIC CLOCKS
We consider a system of mutually delay-coupled electronic clocks, so called phase-locked loops (PLLs) [24,25]. Each PLL of the clock networks considered here consists of a phase-detector (PD), a loop-filter (LF), a voltage controlled oscillator (VCO), an inverter (INV) and a feedback-delay component in the feedback path, see Fig. 1. Heterogeneity in these components manifests in heterogeneous intrinsic frequencies, filter cut-off frequencies, transmission-and feedback-delays. The PLLs are mutually connected, i.e., each VCO sends its output signal x k (t) to at least one other and receives time-delayed signals x l (t − τ kl ) from at least another oscillator in the network. The PD generates a signal that contains information about the phase relations between the input signals x l (t − τ kl ) from other PLLs and the internal feedback signals x k (t − τ f kl ) x PD k (t) = 1 n k N l=1 c kl h PD φ l (t − τ kl ), φ k (t − τ f kl ) , (1) where h PD ( · ) denotes 2π-periodic coupling functions that depend on the type of phase detection at the PD (e.g. XOR or analog multiplier), φ l (t − τ kl ) denotes the phases of the input signals from other PLLs, delayed by transmission-delays τ kl , φ k (t − τ f kl ) denotes the phase of the feedback signal, delayed by a feedback-delay τ f kl , and n k = l c kl is the node degree and denotes the number of the nodes' input signals. The c kl , equal to one if there is a connection between PLL l and PLL k and zero otherwise, denote the elements of the adjacency matrix that specifies the coupling topology. The LF then processes the PD signal which yields the filtered control signal x C k (t) where p k (u) denotes the impulse response of the LF. Usually the LF is implemented as a RC low-pass. Its transfer-function in Laplace-domain is represented in time-domain by the Gamma-distribution [26]. A large class of LFs can hence be modeled as p k (u) = p(u; a k , b k ) [26], where p(u; a, b) = u (a−1) e −u/b b a Γ(a) , and The LF is characterized by the order a k of the filter, and the scale parameter b k which are related to the cut-off frequency of the LF as ω c k = 2πf c k = 1/(a k b k ). The VCO is operated such that it responds linearly to the control signal x C k (t)φ where k = 1, . . . , N indexes the PLLs in the network, ω 0 k denotes the intrinsic frequency and K VCO k denotes the input sensitivity of VCO k. Using Eq. (1), Eq. (2) and Eq. (4) we obtain the phase model . We here consider the high frequency components to be ideally filtered by the LF [13]. In the case of digital signals and an XOR phase detector, as will be shown in the section with the experimental results, the coupling function h[ · ] is a triangular function [5]. In the following, we consider a cosine coupling function associated with analog signals and a multiplier PD, see [13]. Note that for an LF impulse response peaked at zero, p(u) = δ(u), a sinusoidal coupling function h( · ), and τ kl = τ f kl = 0, Eq. (5) reduces to a Kuramoto model of coupled phase oscillators with heterogeneous intrinsic frequencies [10,28].

III. SYNCHRONIZED SOLUTIONS AND LINEAR STABILITY
We now consider a system of two delay-coupled phaselocked loops. This minimal system is suitable to exemplify the implications of component heterogeneities for the dynamics of self-organized synchronization. The general case of N delay-coupled heterogeneous PLLs is shown in the Supplementary material V. The instantaneous frequencies of the two analog delay-coupled PLLs are given bẏ We are interested in phase-locked synchronized states, i.e., all oscillators evolve with same collective frequency Ω and have constant phase-lag β between them φ 1 = Ωt; φ 2 = Ωt + β.
In this case synchronized solutions with global frequency Ω and phase difference β are given by the following transcendental equations (derivation see Appendix A) Ω =ω +K cos(Ωτ e ) cos(B) − ∆K 2 sin(Ωτ e ) sin(B), where B = (Ω∆τ e )/2 + β and The heterogeneous parameters, ω 1,2 , τ 12,21 , τ f 12,21 , K 1,2 , ω c1,c2 are written in terms of their meanx = (x 1 + x 2 )/2 and difference ∆x = x 2 − x 1 . Parameter symbols without a subscript indicate identical parameters for both oscillators. We also definedτ e = (τ −τ f ), and ∆τ e = (∆τ − ∆τ f ). For the case of identical oscillators, in-and antiphase synchronized states exist, characterized by a common global frequency Ω and a constant phase difference β which is either zero or π, respectively. How heterogeneous parameters affect the existence of synchronized solutions, the global frequencies and phase configurations given by Eq. (8) will be discussed in detail in Sec. IV.
The stability of such solutions can be determined by analyzing the response to small perturbations. More precisely, we can look at the perturbation mode related to the phase-difference between the oscillations of two PLL clocks. The characteristic equation governing the exponential growth or decay rate λ of these perturbations is given by (derivation see appendix B) where the λ denote the eigenvalues of the perturbation modes,p(λ) = (1 + λb k ) −1 denotes the Laplace transform of the impulse response function introduced in Eq. (3), and with h being the derivative of the coupling function with respect to its argument. The eigenvalue (λ max = σ + iγ) with the largest real part σ dominates the long-term behavior and determines the stability of synchronized states. For σ > 0 perturbations to synchronized states grow and the solution is unstable, whereas σ < 0 implies linear stability and that perturbations decay at a characteristic time-scale t c = −σ −1 . The case of σ = 0 denotes marginal stability, i.e., perturbations neither grow or decay. This is always a solution to Eq. (10) and represents an equal shift of the phase of each oscillator. The imaginary part γ = Im(λ max ) denotes the frequency of the perturbation response. Since the stability depends on the frequency and phase-configurations, the delays and cut-off frequencies, it can also be modified when heterogeneities are introduced.

IV. EFFECTS OF HETEROGENEITIES ON SYNCHRONIZATION
In this section, we systematically analyze the effects of heterogeneity in the different clock components. The results are compared to those obtained from an experimental setup of two identical delay-coupled PLLs. We show how the notion of in-and antiphase synchronized states, observed for identical oscillators, becomes blurred when heterogeneities in the intrinsic frequencies, coupling strength, transmission and feedback-delays are introduced. For those cases we define so called asymptoticinphase and -antiphase synchronized states that approach the in-and antiphase synchronized states for identical clocks as the heterogeneities approach zero. In the following sections we will introduce the heterogeneities gradually, starting with heterogeneity in the intrinsic frequencies, followed by the heterogeneities in the other parameters.

A. Heterogeneous intrinsic frequencies
Here we analyze how synchronization in a system of two delay-coupled PLLs is affected by heterogeneous intrinsic frequencies. The results that we obtain for the second-order Kuramoto-model are identical to the ones obtained for the first-order Kuramoto-model with timedelayed coupling and no signal filtering [10]. We assume that signal transmission-delays, coupling strengths and impulse response functions of the LFs for the two oscillators are equal, i.e., τ 12 = τ 21 = τ , K 1 = K 2 = K, p 1 (u) = p 2 (u) = p(u). The feedback-delays are set to zero, τ f 12 = τ f 21 = 0. Under these assumptions we find for the frequencies in Eq. (8) and the corresponding phase differences β FIG. 2: (Color online) Global frequency Ω, the phasedifference β, the perturbation response rate σ and the corresponding modulation frequency γ as a function of the transmission-delay τ for two mutually delay-coupled PLLs with K = 0.25 radHz and ωc = 0.25 ×ω radHz. The left column shows the results with identical (ω1,2 = 2π radHz) and the right column the results for heterogeneous intrinsic frequencies (ω1,2 = (1 ∓ 0.02)2π radHz). The blue (dark gray) and red (light gray) curves correspond to the inphase and antiphase (∆ω = 0) or asymptotic-inphase and -antiphase (∆ω = 0.04×2π radHz) synchronized states. The thick curves denote stable solutions (from Eq. (14), with σ = Re(λmax) < 0) and the thin curves denote unstable solutions.
This quantifies the phase differences β as a function of the detuning ∆ω between the intrinsic frequencies. The phase differences are no longer 0 (inphase) or π (antiphase) as in the case of identical oscillators. Instead, new phase configurations emerge due to the detuning. They depend on the ratio of frequency differences to coupling strength ∝ (∆ω/K), see Eq. (12). The detuning also implies, see Eq. (11), that synchronized solutions exist only if Therefore, synchronized solutions do not exist below a critical value of the coupling strength K c = ∆ω/2 for any delay value, that is, if |∆ω/2K| > 1. From Eq. (10) the linear stability of such synchronized solutions is given by where α 12,21 = −K sin(−Ωτ ± β) andp(λ) = (1 + λb) −1 . In Fig. 2 synchronized solutions in a system without heterogeneity (left column) and with heterogeneity (right column) in the intrinsic frequencies are shown. The coloring identifies asymptotic-inphase and -antiphase states, blue and red, respectively. In the right column of Fig  for heterogeneous frequencies, we observe regimes for which no synchronized solutions exist, since Eq. (13) is not satisfied. These gaps (windows) appear due to the detuning and are absent when ∆ω = 0, see left column.
In that case, we find that pairs of synchronized solutions go through a saddle-node bifurcation as the magnitude of the detuning ∆ω is increased from ∆ω = 0, see Supplementary material Fig. 20. Furthermore, for identical oscillators, non-generic solutions with coinciding frequencies exist for specific parameter values Ωτ e = nπ, n ∈ N + 0 , see Appendix A. These solutions split up due to the symmetry breaking, i.e., additional solutions emerge as the intrinsic frequencies are detuned. The range of possible global frequencies decreases. For identical oscillators, the global frequencies lie in the range Ω ∈ [ω − K, ω + K], whereas for detuned frequencies this range decreases by ∆ω, i.e., Ω ∈ [ω −K +∆ω/2,ω +K −∆ω/2]. That is, the system can only self-organize to synchronized states with global frequencies that can be reached by all oscillators in the network, see experimental results in [5].
Linear stability of such synchronized states is also affected by the detuning as it depends on the values Ω and β of the synchronized solutions. In Fig. 3 we provide an overview of existence and linear stability of synchronized solutions in parameter space, i.e., the (K, τ )-plane.

B. Heterogeneity in transmission-delays
Now we add heterogeneous transmission-delays, i.e., τ 12 = τ 21 to the system with detuned intrinsic frequencies. The frequencies and phase configurations then are and see Eq. (8). As in the last section we set feedbackdelays to zero, and study the effects of heterogeneous transmission-delays. The stability is given by the characteristic Eq. (10) where α 12 = −K sin(−Ωτ 12 + β), α 21 = −K sin(−Ωτ 21 − β) andp(λ) = (1 + λb) −1 . These equations show that the frequencies of synchronized states depend only on the mean delay valueτ = (τ 12 + τ 21 )/2 but not the delay difference ∆τ = τ 21 −τ 12 . The same is true for the characteristic equation which only depends on the mean delayτ , see [29]. Therefore the global frequencies of synchronized states and their stability remains the same as long as the mean delay valueτ is unchanged, see Fig. 4. However, the effect of heterogeneous delays becomes significant for the phase-differences β with its linear dependence on ∆τ . As a result, the phase difference changes monotonically with ∆τ without affecting the frequency and stability of synchronized states. This is an interesting results since it allows the phase difference to be tuned to arbitrary values at a fixed frequency. In Sec. V we show this in experimental results for delay-coupled digital PLLs.

C. Heterogeneity in signal filtering parameters
Analyzing Eq. (8), we find that the filtering process has no effect on the frequencies and phase differences of synchronized solutions. However, the stability of synchronized solutions depends on the cut-off frequencies ω c 1,2 and the order a 1,2 of the LFs via the expressionp i (λ) = (1+λ/(a i ω c i )) −ai for i = {1, 2}, see Eq. (10). We consider filters of first order (a i = 1) and set ω c 1,2 =ω c ∓ ∆ω c /2. Note that for zeroth order filters the phase equations reduce to a first order Kuramoto model with delayed coupling, see Supplementary material IV A. For delaycoupled identical Kuramoto oscillators, Yeung, Earl and Strogatz [30,31] provide a necessary and sufficient condition for the stability of inphase synchronized solutions, given by Kh [−Ωτ ] > 0. However, in coupled PLL systems, first or higher order LFs have significant effects on stability. In these cases, the stability criterion for systems without a LF is still sufficient but not necessary (Kh [−Ωτ + β] < 0 ⇒ σ = Re(λ) ≥ 0), since additional instabilities arise due to time-scale introduced by the filtering process [13].
In Fig. 5, we illustrate the effects of heterogeneous cut-off frequencies on the stability of one branch of the asymptotic-inphase solutions for different values of K and τ . The mean cut-off frequencyω c = (ω c 1 + ω c 2 )/2 is varied from a minimum value 0.05ω toω. With the meanω c and the difference ∆ω c , the individual cutoff frequencies of the filters are ω c 1,2 =ω c ∓ ∆ω c /2. The magnitude of difference between cut-off frequencies |∆ω c | is varied from zero (identical filters) to a maximum (∆ω c ) max = 2(ω c −0.01ω), such that the lower of the two individual cut-off frequencies (say ω c 1 =ω c − |∆ω c |/2) is not smaller than 0.01ω. The normalized difference ∆ω c /(∆ω c ) max is plotted on the y-axis. We find that the stability depends on the mean and the difference between the cut-off frequencies. For identical LFs (ω c 1 = ω c 2 ) synchronized states are unstable for cut-off frequencies that are small compared to the mean intrinsic frequency, i.e., large integration time of the filters. However, as heterogeneities are introduced to the LFs, i.e., non-identical cut-off frequencies, unstable regimes can become stable. Hence, specific differences ∆ω c in the cut-off frequencies may stabilize synchronized solutions that where unstable in the case of identical LFs. The cut-off frequency induced instabilities also depend strongly on the value of the coupling strength K and the transmission-delay τ . The gray-region indicating unstable solutions becomes considerably smaller as the coupling strength K is decreased. This is interesting as in the case without delays and filtering the stability is usually enhanced with larger coupling strength. The time-scale associated to the loop-filter integration makes the system inert, i.e., effectively responding to outdated information. As the coupling strength increases, we observe that the instabilities introduced by this inert dynamics are enhanced. This is related to the perturbation decay dynamics which is underdamped for most values of the transmission delay, see e.g., values of γ in Fig. 2. The stronger the VCO reacts to the control-signal, the more it tends to overshoot when close to a synchronized state.

D. Heterogeneity in feedback-delay and interplay between feedback, transmission and processing-delays
So far we have discussed the effects of heterogeneities in the intrinsic frequencies of the VCOs, the transmissiondelays and the cut-off frequencies of the LFs. The filtering process is associated with a processing time-delay which is inversely proportional to its cut-off frequency τ c = b = 1/ω c , the LF integration time. As these pro-FIG. 6: (Color online) Global frequency Ω, phase-difference β, perturbation response rate σ and the corresponding frequency γ for two delay-coupled PLLs with heterogeneous intrinsic frequencies ω1,2 = (1∓0.02)×2π radHz as a function of mean feedback-delayτ f with a constant feedback-delay difference ∆τ f = 0.2 s. The blue (dark gray) and red (light gray) curves correspond to asymptotic-inphase and -antiphase synchronized states, respectively. Thick lines denote stable and the thin lines unstable solutions. The results of synchronized solutions with homogeneous feedback-delays ∆τ f = 0 s are shown in gray. Here, transmission-delays are equal and set to τ = 2.0 s, coupling strength is K = 0.50 radHz and the cut-off frequency ωc = 0.25 ×ω radHz.
cessing time-delays increase, perturbation decay rates become smaller [13]. In this section, we discuss how heterogeneous feedback-delays and its interplay with different delay-time scales (transmission and processing-delays) affects synchronized states. The frequencies only depend on the mean of the transmission and feedback-delays (τ , τ f ) while the corresponding phase configurations depend on their differences ∆τ and ∆τ f . The collective effect of these delay-times can be understood in terms of an effective mean delayτ e =τ −τ f and an effective delaydifference ∆τ e = ∆τ −∆τ f . If the feedback-delay is equal to the transmission-delay,τ e will be zero. In that case, the dynamical properties of the system change qualitatively, as the multistability induced by the transmission- In Fig. 6 we plot synchronized solutions and their corresponding eigenvalues as a function of mean feedbackdelayτ f , given the individual feedback-delays as τ f 1,2 = τ f ∓ ∆τ f /2 with fixed difference ∆τ f = 0.2 s. We focus on the effects of feedback-delays and therefore consider equal transmission-delays fixed at τ = 2 s. In this case there is a symmetry of the frequencies aroundτ f = 2 s, i.e.,τ e = 0. As can be seen from Eq. (18) this symmetry aroundτ e is a generic feature as the cosine and sine squared are symmetric around 0. The phase configurations are also modified as in the case of transmissiondelays due to the heterogeneity ∆τ f , however with opposite sign. The characteristic equation (10) for nonzero τ f is For the characteristic equation no single effective delay parameter representing the combined effect of transmission and feedback-delays can be defined. We find that the transmission-delays affect the eigenvalues only through their mean valueτ . Fig. 7 and Fig. 8 show that ∆τ f andτ f have significant effects on the stability of synchronized solutions. As the feedback-delay is increased, synchronized states become unstable. Heterogeneity in the feedback-delays can stabilize and destabilize synchronized solutions that were unstable or stable for identical feedback-delays, respectively.

E. Heterogenous coupling strengths
The coupling strength K denotes the sensitivity of the VCO to the control signal. It interacts with the effects of the transmission-and feedback-delays as well as the filtering process on the dynamics of the coupled system. For example from Ref. [10] it is well-known that the effects of time-delays are more pronounced for larger coupling strengths, see e.g., Fig. 4, and in the Supplementary material Figs. 21-23, for how the number of synchronized states increases with larger coupling The intrinsic frequencies are equal ω1,2 = ω = 1 × 2π radHz and the cut-off frequency is ωc = 0.25 ×ω radHz. The blue (dark gray) and red (medium gray) curves correspond to the asymptotic-inphase and -antiphase synchronized states. The thick curves denote stable solutions and the thin curves denote unstable solutions. The limits of the frequency-range which is accessible to both oscillators with different coupling strength are shown by dashed lines. For comparison, the results for equal coupling strengths, i.e., for ∆K = 0 are plotted with gray (light gray) curves.
K. Simultaneously, more synchronized solutions become unstable as the instabilities introduced by the filtering process become more prominent. On the other hand, the coupling above the critical coupling strength enables synchrony in detuned oscillators by forcing them towards a common frequency, and therefore counters the effects of detuning ∆ω. These effects are modified when PLLs with heterogeneous sensitivity interact with each other, i.e., their coupling strengths are heterogeneous, see Eq. (8). Note that for identical coupling strengths and intrinsic frequencies, there are non-generic solution for specific parameter values when sin(Ωτ e ) = 0. These non-generic pairs of solutions with coinciding frequencies split-up as a result of the symmetry break- ing, i.e., K 1 = K 2 and the system does not become underdetermined. For such heterogeneity in the coupling strengths, the phase-differences and frequencies of the synchronized states are significantly modified, see Fig. 9 and 10. For ∆ω = 0, the maximum deviation |Ω − ω| is now determined by the smallest of the coupling strengths, e.g., Ω ∈ [ω − K 2 , ω + K 2 ] in Fig. 9. For detuned PLLs the range of global frequencies is given by Ω ∈ [max(ω k − K k ), min(ω k + K k )], where k = 1, 2. The stability of these solutions and is given by Eq. (10). The perturbation decay rate depends non-linearly on the difference between the coupling strengths, see Fig. 11. For differences in the coupling strength where one of them approaches zero and the otherK, the decay rates approach the case of perturbation decay in entrained oscillator systems. However, the maximum decay rates are found for values of ∆K that do not equal to the case of entrainment ∆K max or identical coupling, ∆K = 0. For these values where perturbations decay fastest, the decay changes qualitatively from oscillatory to overdamped. This happens as the change in ∆K causes different solutions of the characteristic equation with zero and non-zero imaginary parts to interchange in terms of the maximum real part σ. Fig. 10 suggests that around ∆K = ∓2 the system undergoes a pitchfork-bifurcation and reverse. As soon as detuning of the PLLs is introduced, this changes and something closer to a transcritical bifurcations arises and frequencies of the different synchronized states change gradually with ∆K until they crash onto the other branch. That means that a slight detuning of the intrinsic frequencies of the PLLs can cause an abrupt change in the frequencies of synchronized states.

V. EXPERIMENTS AND SPICE SIMULATIONS
To validate our theoretical predictions, we performed experiments using prototype DPLLs. These are coupled via a microcontroller that buffers and thereby delays their output signals [5]. The DPLL parameters are given in Table I. In the setup currently available we can tune the transmission delay via the microcontroller and the cut-off frequency using the tunable resistance of the RC loopfilter. The coupling-strength and the intrinsic frequencies are fixed system properties whose heterogeneity has been measured and is known [5].  Feedback-delays, inverters and dividers are not available in this prototype setup. A detailed description of how the experiments were carried out can be found in the Supplementary material, section III. We also performed circuit-level Spice simulations using the freeware implementation LTspice [7,32]. In this  framework, the voltages, currents and delay-times associated with the components of the electronic architecture can be simulated, see Fig. 17 in the Supplementary material. The electronic components in the Spice simulation can be configured to match those of the experimental setup by tuning individual voltages, resistances, propagation-delays and capacitances.
Taking into account the measured heterogeneity in the intrinsic frequencies and coupling strengths of the DPLL prototypes we show in theory, experiment and simulation that the phase-model can precisely predict the phasedifferences and global frequencies of synchronized states, see Fig. 13. We also show how the phase-difference depends linearly on the difference in the delays, while the frequency of the synchronized state remains constant as predicted in subsection IV B, see Fig. 14. Furthermore we show an example, how synchrony can be recovered for an unstable synchronized state by tuning the cut-off frequencies to heterogeneous values while keeping the mean cut-off frequency constant, see Fig. 12. This is achieved by tuning the resistance trimmers of two coupled PLLs at runtime from identical values to values that correspond to detuned cut-off frequencies while the mean cut-off frequency is kept constant. As a result the synchronized state becomes stable.

VI. DISCUSSION AND APPLICATION
In this work we studied how breaking parameter symmetry affects the synchronization dynamics in systems of mutually delay-coupled oscillators. We focused our analysis on the tractable case of two oscillators. For studying larger systems we provide the analytic expressions, but their analysis is beyond the scope of this paper. How the effects found for the two PLL system change with, e.g., coupling topologies, the number of oscillators, or in dependence of the distributions of the heterogeneities are interesting starting points for further research and need to be addressed when moving into application. Certainly the control of self-organized synchronization will be highly relevant for the design of large distributed networks of mutually delay-coupled clocks. The selforganization approach to synchronization discussed here is fundamentally and qualitatively different from hierarchical approaches. There is no accumulation of phaseerrors between the oscillators and signaling delays do not necessarily cause phase-offsets between the oscillators in synchronized states, to name two key strengths of self-organized synchronization. However, the complex dependencies of the dynamics on the system parameters pose a challenge when it comes to the control and robustness of synchronization in the light of process, voltage and temperature (pvt)-variations, relevant to modern microelectronics. Contrary to the expected hindering effects of component heterogeneity, known from synchronization in hierarchical networks, we here showed that clock heterogeneities can actually enhance synchronization in networks with flat hierarchy. Specifically, this reflects in the optimal difference of the coupling-strengths that optimizes the decay of perturbations to synchronized states, and heterogeneous LF cut-off frequencies that can stabilize synchronized states. Furthermore time-delays that usually complicate the synchronization of spatially distributed clocks are essential for the stability of synchronized states and control the clocks' mutual phasedifferences. Interestingly this dependence of the phaseconfigurations on the effective delays is linear.
Hence, the complex interplay between the different system parameters leads to rich dynamics and can provide the means to tune into synchronized states with, e.g., arbitrary phase relations. These findings have strong implications for the application of self-organized synchronization in engineering. Such an approach without hierarchical structures can only be taken if synchronization dynamics can be controlled. As heterogeneities in the clock components are inevitable in modern integrated circuitry, adding transmission and feedback-delay elements and allowing for tunable loop-filters can enhance the control of the dynamics and perturbation decay processes in such systems. Defined phase differences between frequency synchronized clocks for example enable beamforming, i.e., spatially focused emission of electro-magnetic waves, relevant for communication, audio and radar applications [34][35][36][37]. A common, spatially distributed time-reference with locked phases can be the enabler of precise indoornavigation systems, as the accuracy of the localization depends directly on the synchrony of the satellite-nodes. Servers for globally available data-bases could be kept in sync to reduce the time-uncertainty of their timestamps [1,2]. Furthermore, cut-off frequencies much smaller than the mean intrinsic frequency make such clock networks inert and lead to slow perturbation de-cay [13]. At the same time however, such small cut-off frequencies will enhance the signal-to-noise ratio within the circuitry [38]. These properties, together with the results that heterogeneity in the cut-off frequency can stabilize synchronized states that would have otherwise been unstable, show the potential for optimization of synchronization in large networks of spatially distributed electronic clocks [39,40]. We envision that the effects studied here will be used in next-generation synchronization layers for large spatially distributed systems. All free system parameters that can be tuned as the system operates would then be used to steer self-organized synchronization dynamics towards application-specific requirements within the same conceptual setting.

VII. ACKNOWLEDGMENTS
We thank Johannes Fritzsche, Frank Jülicher, David Jörg, Rabea Seyboldt, Kevin Bassler, and Lennard Hilbert for stimulating discussions. This work is partly supported by the German Research Foundation (DFG) within the Cluster of Excellence 'Center for Advancing Electronics Dresden'.

Appendix A: Synchronized solutions
In the following, we obtain analytic expressions for global frequencies and phase differences in synchronized states for two delay-coupled PLLs with heterogeneous parameters. We consider cosine coupling function, for which the phase dynamics is given by Eq. (6) as, Here the indices i, j = 1, 2; j = i; we will use this convention for i, j here on unless specified otherwise. We are interested in phase-locked synchronized states, i.e., when the oscillators evolve with same collective frequency Ω and have a constant phase-lag β between them, which leads to the ansatz Using this ansatz in Eq. (A1) one has Rewriting this equation in terms of mean-couplingK and coupling-difference ∆K, this equation reads, Therefore, the phase difference β in the synchronized state is given by, We rewrite this expression as  One can obtain the frequencies of synchronized state from above equation by substituting the value of B from Eq. (A7). In sec. IV, equal coupling strengths are assumed to study the effects of heterogeneous intrinsic frequencies, filters and different delay-times. For this case (K 1 = K 2 = K), using Eq. (A3), the expressions for synchronized solutions can be obtained as following.
Multiplying by (τ − τ f ) and using Ω(τ − τ f ) = nπ, we have Due to symmetry in the coupling function, above criteria is fulfilled for two distinct values of β with coinciding global frequency. This degeneracy is lost and the solutions split up when symmetry is broken due to parameters heterogeneity. For generic solutions, we require sin Ω(τ −τ f ) = 0 and obtain the phase differences β that are allowed in the synchronized states. From Eq. (A13) the phase difference β is A transcendental equation for the collective frequency as a function of the phase shift β, the mean delaysτ ,τ f , the mean intrinsic frequencyω and the detuning ∆ω can then be obtained from Eq. (A15) and Eq. (A14) (A16) The generalization of this analysis to finite systems with N oscillators and arbitrary 2π-periodic coupling functions is provided in the Supplementary material, see V.

Appendix B: Stability of the synchronized solutions
In this section we study the linear stability of the self-organized synchronized states characterized in the previous section. These solutions to Eq. (A1) are characterized by common collective frequencies and constant phase differences between the oscillators. Modifying Eq. (A2), we ask how the system responds if subject to small perturbations q 1,2 (t) φ 1 (t) = Ωt + q 1 (t), φ 2 (t) = Ωt + β + q 2 (t). (B1) With this ansatz and using Eq. (A1) we obtain the dynamics of weakly perturbed synchronized states Here we introduced the notations q Expanding the coupling function to first order, the perturbation dynamics can be separateḋ Then, using the exponential ansatz, i.e., q i (t) = d i exp(λt) we obtain where and This expression can be written in matrix form as where, ζ = 1, vector d = (d 1 , d 2 ) T and the elements of the matrix G 2×2 are for k, l = 1, 2. From Eq. (B6), we identify ζ (= 1) as the eigenvalue of G which must satisfy the expression From that we obtain the condition which leads the following characteristic equation for the complex eigenvalues λ = x + iy From this equation we calculate the complex eigenvalues λ = x + iy and determine stability of the synchronized solutions characterized by collective frequency Ω and phase difference β, as discussed in section III.

Supplementary material
Supplementary Material

A. Experimental results from prototype setup
The experimental results were obtained with a prototype setup of mutually delay-coupled digital phase-locked loops (DPLLs). These DPLLs are made of CD4046B integrated circuit-elements [1] with XOR phase detectors and a first-order RC loop-filter with a tunable Vishay Model 63 cermet trimmer resistor ranging from R min = 100 Ohm to R max = 2 MOhm with a tolerance of ±10%, see data-sheet [2]. The transmission delays are controlled using a Digilent ChipKit Max32 microcontroller [3]. Since generating significant transmission-delays in the milliseconds regime for the system of coupled PLLs with operational frequencies in the kilohertz regime would require connection lengths of the order of hundreds of kilometers, we use the microcontroller to buffer the signals of the PLLs and delay their output as specified. The DPLL signals, phase-differences, transmission-delays and intrinsic and global frequencies were then measured using a PicoScope 2205 mixed-signal oscilloscope [4] and via the buffer history of the ChipKit Max32 microcontroller. The same experimental setup had been used in [5].
The experimental results presented in, e.g., Fig. 13 and 14 were obtained as follows. In a first step we measured the transmission delay-times between the PLLs and obtained the mean and the difference for each bidirectional connection. After transients had decayed, we captured 50 digital waveforms of 10 ms, each with a sampling-rate of 1.22MHz. The frequency Ω j of the output signal can be calculated from the time-period T Ωj , i.e., the time-interval between two consecutive rising (or falling) edges of measurement j, where j = 1, . . . , J T indexes all measurements in all 50 waveforms that were recorded, see Fig. 15. From these measurements we obtain sets of instantaneous Ω j and β j values that are calculated from consecutive rising edges in each waveform for all 50 waveforms. From this data-set, we calculate their respective mean values Ω, β and the corresponding errors. In order to calculate an average phase-difference β, we use the complex order parameter defined as r β exp(iβ) = (J T ) −1 j exp(iβ j ). As for large J T the value r β approaches zero if the angles β j are equally distributed between [0, 2π) and one if all β j 's are equal, the value (1 − r β ) was used to calculate the measurement error of β. The error can then be specified in the range [0, π) by multiplying the value (1 − r β ) by π. Note that the incomplete cycles at the beginning and at the end of the signals and the transients after the system had been turned on were excluded.

Heterogeneity in intrinsic frequencies and coupling strengths
Using the results of the analysis of the phase model for mutually delay-coupled heterogeneous DPLL elements, we study the global frequencies and phase-differences as a function of the transmission delay. The time-series of the signals were recorded and measured in a system of two mutually delay-coupled DPLL prototype units with detuned intrinsic frequencies and coupling strengths, see also table I in section V of the paper. From these measurements for different transmission-delays we obtained the times of rising and falling edges which we use to extract a phase.
Here we made the assumption that the frequency is constant between two consecutive rising or falling edges and advance the phase linearly by 2π in the resulting time window, see Fig. 15. Using this phase time-series we can then calculate the instantaneous frequencies of the DPLL elements, the global frequency of synchronized states, the Kuramoto order-parameter and the phase-differences. The relation between the perturbation decay-rates, frequencies and the magnitude of the order-parameter R(t) is given by where σ = Re(λ), γ = Im(λ) and C 0 a constant related to the initial perturbation vector [6]. Hence, fitting the time-evolution of the Kuramoto order-parameter allows to measure the perturbation decay-rate and in the case of underdamped dynamics its associated frequency, which can then be compared to the theoretically predicted values. However, due to the limitations of the current experimental setup the perturbation decay-rates could not measured here and will be addressed in the next-generation experimental setups that are planned. The experimentally obtained results with the heterogeneous elements are in very good agreement with the values predicted by the phase model, see Fig. 13 in section V of the paper.

Heterogeneity in transmission-delays
The transmission-delays are controlled using the microcontroller. Using the internal clock of the microcontroller the output signals are delayed by defined times and then sent out to the DPLLs in the network. This can be done individually for each transmission-delay, i.e., for the signal transmission-delay from PLL 1 to PLL 2 and vice versa. For the experiments, the delays are changed such that the mean value of the transmission-delays remains constant. For each value of delay-difference we measure the individual transmission-delays, record the signal timeseries and the frequencies. To confirm the value of transmission-delay set by the microcontroller we measured the time-difference between the same rising edge at the output of each PLL and the input at its coupling partner using the oscilloscope. We plot the theoretical predicted values for the global frequencies and phase differences for DPLL elements with heterogeneous intrinsic frequencies and coupling strength. The experimental results and corresponding LTspice simulations are in very good agreement with the results predicted by the phase model, see Figs. 13-14 in section V of the paper.

Heterogeneity in cut-off frequencies
The current experimental setup allows only a coarse-grained setting of the cut-off frequencies of the loop filter elements using tunable resistors. These potentiometers range from R min = 100 Ohm to R max = 2 MOhm with a tolerance of ±10% and can be tuned using a screwdriver, thereby changing the cut-off frequency of the loop filter. The cut-off frequency is given by where the conductance of the RC loop-filter in our experimental setup is fixed at C LF = 22nF and R LF denotes the resistance. We then measured signal time-series for the case of identical and different cut-off frequencies of the loop- filters, while we kept the mean cut-off frequency constant. Transient measurements were also carried out, detuning the cut-off frequencies of the loop filters while the system was running. In Fig. 12, we show the results of these experiments, i.e., how heterogeneity in the cut-off frequencies of the loopfilters can be used to recover synchronization. For the case presented, the transmission delays τ 12,21 were fixed at 0.7512 ms. We changed the cut-off frequencies ω c 1,2 = (R LF C LF ) −1 by adjusting resistance R LF of the potentiometer. In this setup, we collected the digital waveforms of the two PLLs for 10 s with a sampling-rate of 1MHz and performed 20 such measurements. Initially, the resistances were set to R 1 LF = 131.2 kΩ, R 2 LF = 131.0 kΩ, i.e., loop-filters had approximately equal cut-off frequencies (ω c 1 , ω c 2 ) = (0.0551396, 0.0552238) × 2π radkHz. At the beginning of our measurements, the DPLLs were desynchronized and the synchronized state was unstable. After approximately 4 s, we changed the resistances of the potentiometers to R 1 LF = 490.0 kΩ and R 2 LF = 75.6 kΩ such that the the new cut-off frequencies were ω c 1 = 0.0148 × 2π radkHz and ω c 2 = 0.0957 × 2π radkHz respectively. This change increased the heterogeneity in the cut-off frequencies while the mean value remained approximately the same (ω c = 0.05524 × 2π). As a result, the system synchronized to a state with a common global frequency, see Fig. 12 in section V of the paper.

B. LTspice simulations
We validate the theoretical results obtained from the phase model using LTspice simulations [7]. This uses circuit level models for the dynamics of voltages, currents and processing times of standard electronic components, most of which are available on the market [7]. The operational amplifiers used in the adders are ADA4857 components [8]. As a VCO we use the LTC6900 oscillator [9]. For this device we obtain the input sensitivity and measure the operation frequency as a function of the input voltage V bias , see Fig. 16. From this we then obtained the input sensitivity of the VCO as K VCO = 0.088291MHz/V. The frequency of the VCO depends of the offset voltage V bias as Then we use a setup of two mutually delay-coupled DPLLs to simulate self-organized synchronization on the circuit level, see Fig. 17. We measure and monitor the time-series of the VCO's output voltages, the control signals, the loopfilter signals and the perturbation signals. The perturbation signals are added to the control signal of one of the DPLLs and cause phase-deviations from the synchronized state. Using the time-series of the output-signals of the VCOs we can obtain the phases, phase-differences and the instantaneous frequencies. From these quantities the Kuramoto order-parameter can be calculated as a measure of the synchrony in the system. Also, the perturbation decay rates can be measured from the time-series of the phase-differences or the time-evolution of the order parameter [6]. A network of two DPLLs is characterized by the intrinsic frequencies ω k and the coupling strengths K k of the VCOs, the transmission-and feedback-delays τ kl , τ f kl , and the cut-off frequencies of the loop filters ω c k . In the LTspice model, these parameters are set by different electronic components. The cut-off frequencies are set using the resistors R LF, k of the RC-circuit. The coupling strengths are set by gains in the adder of the XOR signals (S k1 ) and the adder for the bias voltage to the loop-filter signal (S k2 ). These gains S k are set by the resistors R MXOR, k and R PLL, k , where we use the first gain to set the coupling strength and the second gain compensates for damping in the low-pass filter. We derived the corresponding relations of the system parameters and the LTspice model and found for the total and individual gains where A denotes the output amplitude V high of the XORs, and k = 1, 2 indexes the DPLLs. The resistance-values of the resistors of the loop-filter cut-off frequency ω c k and gain S k1 of the phase-detector signal can be calculated from where C LF is the capacitance of the loop-filter capacitor.

IV. NUMERICAL SIMULATION OF THE PHASE MODEL
A. An ath order filter: PLL dynamics as (a + 1)th order differential equation In the following, we show that system dynamics with ath order filter can be written as (a + 1)th-order differential equation. For this, consider the dynamics of N delay-coupled PLLs, rewriting this equation asφ where, (p * x PD k ) represents the convolution operation on p and x PD k . Taking the Laplace transform (L) of this equation, and using the property L [(f * g)(t)] = L(f (t)) · L(g(t)), we get, Here hats (ˆ) represent the Laplace transforms of the respective functions. Sinceφ k (s) = sφ k (s) − φ k (0), one can rewrite this equation as, Laplace transform of the impulse response functionp(s) for ath order filter is given by, Substituting this value in Eq. (4), we obtain, Therefore, with ath order filter, dynamics of the system in Laplace space is given by Substituting binomial expansion for (sb + 1) a i.e., we get, Using the property that the Laplace transform of nth-order derivative of a function (say f [n] (t)) is Eq. (10) reads, Inverse Laplace transform of above equation back into the time domain gives, where, Here δ(t) represents the Dirac-delta function and δ [n] (t) its nth derivative. Hence, the dynamics of the system with ath order filter is governed by (a+1)th order differential equation, which is given by Eq. (13). Few examples (a = 0, 1, 2, 3) are given below. Filter of order zero (a = 0): For a zero-th order filter, P δ k (0) = 0, and the frequencies of the VCOs are given by a simple first order differential equation, Filter of order one (a = 1): Now for a first order filter i.e., a = 1, [1] (t) + φ [1] k (0)δ(t) + b φ k (0) · δ [1] (t) + ω k · δ(t) Therefore, the phase dynamics for the first order filter is given by Filter of order two (a = 2): For second order filter [2] (t) + φ [1] k (0) · δ [1] (t) + φ [2] k (0) · δ(t) [2] (t) + ω k · δ [1] (t) , The dynamics is govern by the third order differential equation Therefore, depending on the order of the loop-filter a, one can rewrite the integro delay-differential equation for the phase-dynamics of coupled PLLs system as (a + 1)th order delay differential equation.
B. Simulation results with second order differential equation For our study we consider the PLLs with first order loop filter. For this system, the phase dynamics can be written as see Eq. (16). The phase detector signals x PD k (t) for analog and digital PLLs are, x PD,analog x PD,digital Therefore, the dynamics of delay-coupled analog PLLs system is given by We use this second order differential equation for our simulations and examine different dynamical scenarios exhibited by the system at different parameter values. One can choose the parameter values and the history for initial conditions and see how the system settles into a modified-inphase, modified-anti-phase states, or remains desynchronized (either synchronized states are unstable or it do not exist). Examples of these cases are shown in the Supplementary Figs. 18 (a)(b)(c) and (d) respectively. These simulations are helpful in quickly verifying the analytical results and to get insight into the behavior of the system. The results shown in Supplementary Fig. 3 indicate the possibility of different dynamical scenarios in parameter space. Here, initial conditions play important role due to the multistability. This behavior can be examined numerically as in Supplementary Fig. 19, where we show the simulation results for the asymptotic phase dynamics in the (K, τ ) plane. In the multistable regions where both modified in-and antiphase solutions coexist, simulations show how the final state of the system depends on the initial conditions. As multiple stable synchronized states may coexist in the system, simulations are also helpful to examine the basin of attraction for these multiple solutions.

V. EXTENSION TO THE NETWORK OF N DELAY-COUPLED OSCILLATORS
Sec. III in the paper provides the analysis for two delay-coupled oscillators with heterogeneous components, where we obtain the expressions to evaluate synchronized solutions and their stability. These transcendental equations can be solved numerically to calculate the global frequency, the phase difference and the corresponding perturbation response rates. For a larger network of delay-coupled PLLs (N > 2) with heterogeneous components, the analysis can be extended as presented in the following.

A. Synchronized solutions
The dynamics of N delay-coupled PLLs, in a general form, is given by, Synchronized solutions are given by the ansatz where Ω is the common collective frequency and β k is the phase deviation. Putting these solutions in Eq. (1), we get were β kl = (β k − β l ). Eqs. (3) provide condition for the synchronized solutions described in Eq. (2). Note that, there are N equations, see Eqs.
(3), to be solved simultaneously for (N + 1) values, {Ω, β 1 , β 2 · · · , β N }. This can be reduced to N equations simply be assuming β 1 = 0 as a reference. The β-values then represent the phase differences of individual oscillators with respect to the reference oscillator. These N values (i.e., the solutions {Ω, β 2 · · · , β N }) can then be calculated from the following N number of equations.
Therefore, the phase-locked solutions satisfy, A compact form of the condition for synchronized solutions can be obtained by averaging over all oscillators, These equations provide synchronized solutions of a notwork of delay-coupled oscillators for a given coupling function and connection topology. For larger N > 2 it is nontrivial to obtain explicit analytic expressions for synchronized frequency and phase differences. One can use, however, for example, linear approximations for coupling function or numerical methods to estimate these solutions.
B. Stability of the synchronized solutions For a system of N oscillators, we here obtain a general expression for the stability of synchronized solutions. In order to calculate linear stability, we apply small perturbations q k to the synchronized solutions, i.e., φ k = Ωt + β k + q k ; k = 1, 2, · · · N ; (7) and examine the behavior of these perturbations. Putting these values in Eq. (1), we obtain, Expanding function 'h' to the first order, we have Therefore, the perturbation dynamics is given bẏ Assuming exponential variations in the perturbations, i.e., q k = d k e λt , we get For notation simplicity, denoting,c kl = c kl /n k and above equation read, Now usingp k (λ) = ∞ 0 du p k (u)e −λu , we have Therefore we find where, In a matrix form we can write, where, ζ = 1, vector d = (d 1 , d 2 , d 3 , · · · , d N ) T and the elements of the matrix G N ×N are, From Eq. (17), ζ (= 1) is the eigenvalue of the matrix G hence must satisfy the characteristic equation det (G − ζ · I) = 0, Using above condition, for the network of N delay-coupled PLLs, one can evaluate eigenvalues λ of the perturbation dynamics and determine stability of the synchronized states.

VI. SUPPLEMENTARY RESULTS WITH HETEROGENEOUS PARAMETERS
A. The effects of detuning As discussed in the paper, one of the effects of the detuned intrinsic frequencies is the appearance of the region in parameter space where synchronized solutions do not exist, see Fig. 3. Here we examine how synchronized states in the system vanish as the detuning is increased. We find that the pairs of synchronized solutions go through a saddle node bifurcation as the magnitude of the detuning ∆ω is increased from ∆ω = 0. This is shown in Supplementary  Fig. 20, where pair(s) of stable and unstable solution collide at the bifurcation point and disappear. Depending on the delay-values, this bifurcation can occur in one (see first and second column), or more pairs of solutions (third column), leading to the regimes in parameter space where synchronized solutions do not exist.
In Supplementary Fig. 21, the behavior of the synchronized solutions is shown as a function of the transmission delay at a higher value of coupling strength as in Fig. 14 in section IV A in the paper. This figure shows that with increasing coupling strength, the transmission delay induced multistability [10] and instabilities caused by the filtering process become more pronounced. With heterogeneous intrinsic frequencies ∆ω = 0 (right column), we can also observed the pair of additional solutions arising due to the symmetry breaking (see main text), which were absent when ω 1 = ω 2 (left column).

B. Heterogeneous transmission delays
In the paper we discussed how the synchronized states are effected by the heterogeneity in the transmission delays. Specifically, for a fixed mean delayτ , we showed that the global frequencies of the synchronized state remain unaffected while the phase-difference changes linearly as a function of delay-heterogeneity ∆τ . Here we present additional results when the delay-heterogeneity ∆τ is fixed and show the behavior of the system as a function of mean delayτ . In Fig. 22 we plot the frequency and phase differences against the mean delayτ for identical (left column ∆ω = 0) and detuned (right column ∆ω = 0.04 × 2π) intrinsic frequencies of the PLLs with a constant delay difference ∆τ = −0.20 s. We find that the frequencies of synchronized states of the system with homogeneous (τ 12 = τ 21 =τ ) and heterogeneous delays τ 1,2 =τ ± ∆τ remain same for equal mean valuesτ . However, the phase-differences deviate from those in the case of homogeneous delay (gray curves for ∆τ = 0) by −Ω∆τ /2, see Eq. (16) in subsection IV B of the paper.   21: (Color online) Global frequency Ω, phase configuration β the perturbation response rate σ and the corresponding modulation frequency γ of the synchronized solutions plotted as a function of the transmission-delay τ for a system of two mutually delay-coupled PLLs at a coupling strength K = 2.0 radHz. All other parameters are the same as in Fig. 2 in the paper. Additional synchronized solutions are visible for detuned intrinsic frequencies, see right column for ∆ω = 0.04 × 2π radHz.

C. Heterogeneous coupling strengths
In section IV E, we have discussed how the difference in coupling strengths ∆K affects the frequencies of synchronized states and phase difference in identical and detuned PLLs. We also observed how the coupling strength heterogeneity changes the perturbation decay rates of the synchronized solutions while the mean of the coupling values is kept constant (see Fig. 11). In Supp. Fig. 23, we present similar results for a larger value of the transmission delay, Global frequency Ω and phase-configuration β of synchronized states plotted as a function of the mean delayτ with a constant delay-difference ∆τ = −0.20 s for a system of two delay-coupled PLLs with identical ω1,2 = 1×2π radHz (left column) and heterogeneous intrinsic frequencies ω1,2 = (1 ∓ 0.02) × 2π radHz (right column). The blue (dark gray) and red (light gray) curves correspond to the inphase and antiphase (∆ω = 0 radHz) or modified-inphase and modified-antiphase (∆ω = 0.04 × 2π radHz) synchronized states. The thick and thin curves denote stable and unstable solutions respectively. The phase configurations in the synchronized states without transmission-delay heterogeneity (∆τ = 0) are plotted with gray curves for comparison. The coupling strength is K = 0.25 radHz and the cut-off frequency is ωc = 0.25 ×ω radHz.
where many synchronized states coexist. For this case, we find for some of the synchronized solutions, the coupling heterogeneity may change the perturbation decay rates to the extant that it crosses zero value, i.e., the stability of the solutions changes. This means that unstable synchronized solutions present for homogeneous coupling strength can be made stable as the system is tuned to heterogeneous coupling strength.