Self-organized synchronization of digital phase-locked loops with delayed coupling in theory and experiment

Self-organized synchronization occurs in a variety of natural and technical systems but has so far only attracted limited attention as an engineering principle. In distributed electronic systems, such as antenna arrays and multi-core processors, a common time reference is key to coordinate signal transmission and processing. Here we show how the self-organized synchronization of mutually coupled digital phase-locked loops (DPLLs) can provide robust clocking in large-scale systems. We develop a nonlinear phase description of individual and coupled DPLLs that takes into account filter impulse responses and delayed signal transmission. Our phase model permits analytical expressions for the collective frequencies of synchronized states, the analysis of stability properties and the time scale of synchronization. In particular, we find that signal filtering introduces stability transitions that are not found in systems without filtering. To test our theoretical predictions, we designed and carried out experiments using networks of off-the-shelf DPLL integrated circuitry. We show that the phase model can quantitatively predict the existence, frequency, and stability of synchronized states. Our results demonstrate that mutually delay-coupled DPLLs can provide robust and self-organized synchronous clocking in electronic systems.


Introduction
In spatially extended electronic systems that require coordination of a large number of components, it is a challenge to provide a synchronized time reference [1][2][3]. Such systems include, e.g., global positioning systems, servers on the internet, multi-core and multi-processor computer architectures, network-on-chips, and large antenna and sensor arrays [4][5][6][7]. A general issue for synchronization of oscillators are signal transmission delays caused by finite signal propagation speed [8][9][10]. A processor running at 1GHz has a period of 1ns, the time in which a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 signals send at the speed of light can travel 30cm. This implies significant transmission delays at length scales in the centimetre range for such frequencies. In most cases, the time reference in electronic systems relies on hierarchical concepts, i.e., a high quality master-clock feeds forward a synchronization signal into a clock-tree, entraining all other slave-clocks [11]. In such systems, signal transmission delays have to be compensated, e.g., all connections in a clocktree have to be of equal length [11,12]. Such strategies become space and energy inefficient for large systems [13]. Moreover, hierarchical entrainment is vulnerable to fluctuations and failures, e.g., due to noise and cross-talk [11,14]. Hence it is a challenge to achieve precise and robust synchronization of autonomous clocking units [15,16].
Robust self-organized synchronization without hierarchical structures has been found in large collections of oscillators in nature. In biological systems, such as neuronal networks, coupled genetic oscillators, cardiac pacemaker cells, and in flashing fireflies, synchronization is attained in a self-organized way [17,18]. Through mutual coupling of individual oscillatory units, these systems are able to synchronize robustly in highly noisy environments and in the presence of considerable communication delays [19][20][21][22][23][24]. For instance, in neuronal systems, axonal conduction delays are typically in the range of 1-100 of milliseconds [25], and therefore a considerable fraction of the typical neural oscillation frequencies between 1 and 100 Hz. This is achieved in the absence of an entraining master clock. It has been shown that such a decentralized synchronization strategy enables robust self-organization of synchronous states also in analog electronic networks, thereby providing a position-independent time reference [26]. Recently, we have developed a phase oscillator description which captures the dynamics of digital electronic oscillators [27]. We also showed that electronic networks provide excellent experimental conditions to study the nonlinear dynamics of mutually coupled oscillators [28][29][30] in a controlled environment with adjustable system parameters and transmission delays [27].
Here, we study synchronization dynamics in networks of delay-coupled digital phase-locked loops (DPLLs). We show theoretically that such networks operating on digital signals can be formulated in terms of a phase oscillator description and study the existence, stability and the frequencies of different frequency synchronized states, e.g., checkerboard, m-twist and in-phase synchronized states. In contrast to implicit expressions for the frequencies of the synchronized states in analog systems, we here provide closed analytical expressions for the frequencies of synchronized states in digital systems. These can be obtained explicitly due to the piecewise linear coupling function. Furthermore, we discuss the role of the characteristic integration time related to signal filtering within the DPLLs. This integration time leads to memory effects that generate inertia-like behavior. This can lead to transient waves that travel through the network, initiated by local perturbations [31]. Our theoretical results are supported by experiments performed on real networks of off-the-shelf DPLL integrated circuitry. These experiments include networks of up to 9 DPLLs with different coupling topologies, such as nearest-neighbour coupled networks with open or periodic boundary conditions. In Section 1, we develop a phase model of N mutually delay-coupled DPLLs. In Section 2 and Section 3, the existence and stability of different types of frequency-synchronized solutions for networks of mutually coupled DPLLs are investigated. In Section 4, we use our theory to illustrate the behavior of DPLL systems for specific examples. In Section 5, we present experimental data for systems of mutually coupled DPLLs with different network topologies and show that the measurements can be understood quantitatively using our theory. Finally, conclusions are drawn in Section 6.
signal to the feedback of the output signal [32]. Such clocks are widely used in electronics to control and synchronize components with respect to a master-clock signal. The essential building blocks of a DPLL are the phase detector (PD), the loop filter (LF) and the voltage-controlled oscillator (VCO). The PD can be implemented as an XOR gate and compares the external signal x l to the feedback signal x k , see Fig 1. All external signals x l and the feedback signal x k are compared using XOR gates at the input for each input signal individually. High frequency components in the output signal x PD k of the PD are filtered by the LF, a low pass filter allowing only low frequencies to pass. Then the output signal x C k of the LF controls and modulates the VCO's frequency. We consider a DPLL containing an optional signal inverter (INV) between the VCO and the PD (see 1). This inverter enables to control the properties of synchronization by introducing an additional phase shift π, independently of the instantaneous frequencies, as will be shown in Section 5. We first discuss a single DPLL to understand its function. In the following we extend our description to networks of arbitrarily coupled DPLL clocks.

Phase description of a single DPLL
We develop a phase description of a single DPLL receiving a time-delayed signal from a distant source. The components of the DPLL can be represented as follows. The VCO outputs a rectangular oscillation with amplitude 1, where ϕ(t) is the phase and P is the square-wave function with Fourier representation P is 2π-periodic and takes the value −1 for −π < ϕ < 0, the value 1 for 0 < ϕ < π, the value 0 for ϕ = 0 and ϕ = π. The PD compares the external input signal x ext with the signal of the internal feedback x using an XOR operation. We account for transmission delays with a single  delay time τ in the input signal, where È denotes the XOR operation and xðtÞ ¼ 1 À xðtÞ. The output signal x PD of the PD is filtered by the LF whose output signal x C is given by where p(u) denotes the impulse response of the LF and satisfies R 1 0 du pðuÞ ¼ 1. This control signal x C is passed to the VCO and determines the VCO's oscillation frequency according to where ω 0 is the minimal frequency of the VCO for zero input and K VCO is its sensitivity. We here consider a linear frequency response of the VCO, such that the frequency response is proportional to the control signal x C . Fig 2 shows this linear response together with the measured response curve from our experimental setup, discussed in Section 5. We consider the loop filter as an ideal low-pass filter that perfectly damps high frequency components, see A. Using Eqs (3) and (4) in Eq (5), we obtain where, for compactness of notation, we have defined ω = ω 0 + K VCO /2 and K = K VCO /2 and where Δ is the triangle-wave function with Fourier representation Δ is piecewise linear with first derivative Δ 0 (ϕ) = −2/π for −π < ϕ < 0, slope Δ 0 (ϕ) = 2/π for 0

Networks of mutually coupled DPLLs
Networks of N mutually delay-coupled DPLLs can be analyzed using an extended version of the phase model Eq (6). If there is more than one input signal to a DPLL, each input is processed with the feedback by an XOR gate individually. Subsequently all signal paths are joined into an analogue average which yields the PD's output signal x PD (see Fig 3). The phase model for a network of N coupled identical DPLLs and equal transmission delays reads where ϕ k is the phase of oscillator k 2 {1, . . ., N}. The connection topology between all DPLLs is described by the coupling matrix C ¼ ðc kl Þ with c kl 2 {1, 0}, where c kl = 1 indicates a connection between DPLL k and DPLL l. The coupling strength is normalized by the number of input signals, n(k) = ∑ l c kl , generating the average input signal to the LF of DPLL k. Eq (8) describe a network of N coupled DPLLs, taking explicitly into account a filter impulse response p(u) and identical transmission delays τ. Identical transmission delays can, e.g., be achieved in a regular square lattice with nearest-neighbor coupling.
In large systems of delay-coupled oscillators, many different dynamic states can exist, e.g., phase-synchronized and frequency-synchronized states as well as more complex phase patterns such as spirals, traveling waves, and irregular or chaotic states [33,34]. Here, we are interested in synchronized states with a constant collective frequency and constant phase relations between oscillators. In the following two sections, we analyze the existence and stability of two important classes of such states, the in-phase synchronized states (Sec. 2) and the mtwist synchronized states (Sec. 3).

In-phase synchronized states 2.1 Frequency of the in-phase synchronized state
In globally in-phase synchronized states (see Fig 4), all DPLLs evolve with the same collective frequency O and have no phase lag relative to each other, Inserting Eq (9) in Eq (8), we find a condition for the collective frequency, Here we have used R 1 0 du pðuÞ ¼ 1, n(k) = ∑ l c kl , and the fact that Δ is even. An in-phase synchronized state exists if Eq (10) has a solution in O. Explicit solutions to Eq (10) can be obtained graphically as the intersection of the lhs and the rhs of Eq (10), see Fig 5. The intersection points can be expressed explicitly as where the sign ± refers to intersections with the rising and falling slopes of Δ, respectively. The indices j fall within a range j 2 N between A and B − 1/2 for O þ j and between A and B + 1/2 for O À j , where A = (ω − K)τ/2π and B = (ω + K)τ/2π. For a given set of parameters, multiple solutions with different collective frequencies can exist. The collective frequency is independent of the number N of oscillators and the coupling matrix C. The number of solutions depends on  the parameters, in particular on the transmission delay τ. In the following section, we analyze the stability of these states.

Stability of the in-phase synchronized state
In-phase synchronized states are robust against fluctuations if small perturbations decay. Performing a linear stability analysis, we obtain the stability properties of the state given by Eq (9) [35,36]. We consider the dynamics of small perturbations q k (t), defined by The linear dynamics of the perturbation is obtained by expanding Eq (8) to first order with respect to q k , where α = KΔ 0 (−Oτ). Using the definition of the triangle-wave function Δ, Eq (7), α explicitly reads where P is the square-wave function Eq (2). Using the exponential ansatz q k (t) = v k e λt in Eq (13), with complex frequency λ and perturbation v k , we obtain the characteristic equation where the Laplace transform of the impulse response,pðlÞ ¼ R 1 0 du e À lu pðuÞ, is the transfer function of the LF and d kl = c kl /n(k) are the components of the normalized coupling matrix D ¼ ðd kl Þ. We solve Eq (15) for the unknowns λ and v = (v 1 , . . ., v N ). To this end, Eq (15) can be rewritten in the form rðlÞv ¼ Dv, where rðlÞ ¼ e lt ðlpðlÞ À 1 a À 1 þ 1Þ. This reveals that v must be an eigenvector of the matrix D, which satisfies Dv ¼ zv, and where z is the corresponding eigenvalue. Hence, λ must satisfy ρ(λ) = z, which explicitly reads Note that each value of z generates an infinite discrete set Λ z of solutions λ. If D is diagonalizable, arbitrary perturbation vectors can be expressed as linear combinations of the eigenvectors v. The matrix D always has the special eigenvector v = (1, . . ., 1) with z = 1, corresponding to a global phase shift of all oscillators. This eigenvector exists because ∑ l d kl = 1 for all k by definition of D. Since a global phase shift does not perturb the synchrony of the system, we do not consider the corresponding perturbation modes in the following. The long-time evolution of the perturbation is dominated by which is the solution with the largest real part within the set of solutions S z 6 ¼ 1 Λ z [26]. The real part σ of λ 0 is the perturbation response rate. The in-phase synchronized state is stable if and only if σ < 0. The imaginary part β is a frequency that leads to a periodic modulation of the phase. The perturbation q k (t) thus generates side-bands in the oscillator signal with frequencies O ± β [26].

m-twist and checkerboard synchronized states 3.1 Frequency of m-twist and checkerboard synchronized states
We now consider frequency-synchronized states with a fixed phase difference between neighboring DPLLs (see Fig 4). In one dimension, the phase profile of such m-twist states has the form Introducing periodic boundary conditions enforces that the phase difference can only take the values where m 2 {0, . . ., N − 1} is the total number of 2π-phase increments along the ring. The case m = 0 corresponds to the in-phase synchronized state described in Sec. 2. For N even, the case m = N/2 corresponds to a checkerboard pattern, i.e., a phase difference of θ N/2 = π between neighboring oscillators. The in-phase and checkerboard synchronized states also exist in the case of open boundary conditions. Inserting the Ansatz Eq (18) into Eq (8) using the nearest-neighbor coupling matrix in one dimension for the collective frequency, which depends on m. Here we have used

Stability of the m-twist synchronized state
We study the linear stability of the states given by Eq (18) by using the Ansatz where q k is a small perturbation to the m-twist synchronized state. The linearized dynamics of the perturbations q k can be written as where the sum runs over + and −, with α ± = KΔ 0 (−Oτ ± θ m ). Using the definition of the triangle-wave function Δ, Eq (7), α ± explicitly read Inserting the exponential ansatz q k (t) = v k e λt into Eq (22), the characteristic equation reads where we imply periodic boundary conditions. As in Sec. 2.2, this equation can be rearranged and written in vector notation as gðlÞv ¼ Ev, where gðlÞ ¼ e lt ð2lpðlÞ À 1 þ a þ þ a À Þ and the matrix E has entries e kl = α + δ k,l−1 + α − δ k,l+1 + α + δ k,N δ l,1 + α − δ k,1 δ l,N . Hence, for this equation to possess solutions in λ, γ(λ) must be an eigenvalue of E. The eigenvalues η j of E are given by where j = 0, . . ., N − 1. Note that the matrix E and hence the eigenvalues η j depend on both the coupling matrix C and the other parameters through α + and α − . The characteristic equation γ(λ) = η with η being an eigenvalue of E can be written as Eq (26) has a discrete infinite set of solutions for each η.

Stability of the checkerboard synchronized state
So far, we have considered periodic boundary conditions. However, the checkerboard synchronized state, characterized by a phase difference θ m = π between neighboring oscillators, also exists for open boundary conditions. The characteristic equation which determines the linear stability of this state is where α + has been defined in Eq (23) and the eigenvalues η j are given by using the fact that α + = α − for θ = π.

Specific DPLL systems as examples
We now discuss the collective frequency and the stability of the in-phase, 1-twist, and checkerboard synchronized states of three coupled DPLLs. We consider both a ring and a chain coupling topology (see Fig 4).

Collective frequency
The collective frequency O for the in-phase, 1-twist, and checkerboard state, determined by Eqs (10) and (20), depends on the intrinsic frequency ω of the VCOs, the coupling strength K, and the transmission delay τ. Fig 6A shows O for all three states as a function of τ. For a fixed value of τ, multiple states with different collective frequencies can be stable simultaneously. The dependence of the collective frequency on the transmission delay is qualitatively different for the in-phase and the 1-twist synchronized states. The range of possible collective frequencies is larger for the in-phase state.

Stability and perturbation response rates
For σ < 0, the synchronized state is linearly stable and the perturbation response rate −σ is a measure of the time scale of synchronization. This time-scale or equivalently the perturbation response time can be calculated as t c = −1/σ. It corresponds to the time in which perturbations have decayed to q k (t c )/q k (0) = 1/e and can be obtained using q k (t) = q k (0)e λt . Fig 6 shows σ as a function of Oτ for the in-phase and 1-twist synchronized states in a ring of 3 coupled DPLLs (see Fig 6B) and for the in-phase and the checkerboard synchronized states in a chain of 3 coupled DPLLs (see Fig 6C). Note that there are parameter regions for which none of these states is linearly stable, indicating that the system attains more complex dynamical states, whose characterization is outside the scope of the present work. For the special case of no delay, τ = 0,  (16) and (26). Here we consider the large class of loop filters with the transfer function where a is the order of the filter and ω c is the cutoff frequency [37]. Fig 7 shows the perturbation response rate σ for two examples in which the filter transfer functionp is given by Eq (29) with filter order a = 1 and a = 0, where the latter case corresponds to no filtering. Shaded regions indicate stable in-phase synchronized solutions. Note that the regions where the inphase synchronized state is stable differs for these two cases. It has been shown earlier that for the case of no filtering (a = 0), the in-phase synchronized state is stable if and only if α > 0 with α given by Eq (14) [38]. Fig 7 thus shows that in the presence of a filter, linear stability is no longer determined by the sign of α. Hence, the presence of a filter can change the stability of synchronized states. Moreover, filtering also affects the time scales of synchronization as has been shown earlier [26]. Therefore, effects of filters have to be taken into account in models of coupled PLLs in order to obtain correct stability properties. For filters of first order (a = 1), the Self-organized synchronization of delay-coupled DPLLs effect of the filtering can be understood as inertia of the phase variable due to characteristic integration time of the filter, see C.

Experimental measurements on coupled DPLLs
We designed prototypes to conduct experiments with networks of mutually coupled DPLLs to compare the synchronization properties of a real system to our theoretical predictions. The prototype networks were set up with DPLLs (CD4046B [39], specifications given in Table 1 and Fig 8), a Digilent ChipKit Max32 microcontroller [40], and a PicoScope 2205 Mixed-Signal oscilloscope [41] (see Fig 9). We choose the kHz regime for our prototypes to make sampling of the signals simple, and because the effects introduced by the delay only depend on the relation between delay and period of the uncoupled oscillators. If signals can be transmitted with about two thirds of the speed of light [9], clocks operating in the kHz range would require connection lengths of the order of kilometers to achieve significant transmission delays. The necessary transmission delays were achieved by artificially holding back the signals in a microcontroller according to a chosen value of the delay. We measure the output signals of the DPLLs using an oscilloscope. From the data the phase time series of each of the DPLLs of the network was obtained. We use this data to measure the exponential relaxation time to a synchronized state and frequency of the synchronized state. The measurement parameters are given in Table 2. We performed measurements on three different systems, analyzing for each system different synchronized states: (A.) a system of 2 mutually coupled DPLLs, (B.) a ring of 3 mutually coupled DPLLs, and (C.) a two-dimensional square lattice of 3 × 3 mutually coupled DPLLs with periodic and open boundary conditions. Note that system (B.) corresponds to one of the examples discussed in the last section. The specifications of the used DPLLs are given in Table 1. After a system is started, we first let the DPLLs oscillate independently during a randomly chosen time and then turn on the coupling. Because of small differences in intrinsic frequencies, oscillators have an arbitrary phase difference at the moment they are coupled. Subsequently, the system settles into a synchronized state for a wide range of delay times. All experimental results are averages over up to 40 independent measurements with the standard deviation being smaller than the symbol size in all plots presented. The number of data points that are averaged depends on how likely the respective synchronized state is attained starting from arbitrary initial conditions. We compare these measurements to results of our phase model using the average DPLL parameters for the intrinsic frequency ω and the coupling strength K. Self-organized synchronization of delay-coupled DPLLs

Two mutually coupled DPLLs
We first discuss the simplest case of two mutually coupled DPLLs (no. 1 and no. 2 in Table 1). Two types of synchronized states are observed: the in-phase and the checkerboard synchronized state (see Fig 4). Fig 10 shows the measured collective frequency O for different values of the transmission delay, together with the results obtained from the phase model, Eq (10). Which of the two synchronized states is found depends on the value of the transmission delay. For certain ranges of the delay τ both synchronized states coexist and the system is bistable. In the coexistence region, synchronized states are selected stochastically after transient dynamics.
The VCO response becomes nonlinear and saturates in the VCO clipping region (see Fig 2). Outside these VCO clipping regions, where the VCO response is linear, the measured values of the collective frequency are in good agreement with the results of the phase model, see Fig 8. Inside the VCO clipping regions, the linear response approximation of the VCO's dynamic frequency, Eq (5), becomes inaccurate and the measured collective frequency deviates from the theoretical results. Adding a signal inverter in the feedback loop (see Fig 1) exchanges the collective frequency and the stability of the in-phase and the checkerboard synchronized state (compare Fig 10A and 10B). Alternatively, the inverter can also be added to the output of each PLL or each input with the same result. From experimental measurements, we obtained the exponential relaxation time σ −1 of the phase difference between the two DPLLs by fitting an exponential function to the data.  Table 1. Red dots indicate the intrinsic frequency ω, bars indicate the frequency range where the VCO response is linear (see Fig 2). Shaded areas indicate the system's clipping regions; outside these regions the response of all VCOs is linear.
doi:10.1371/journal.pone.0171590.g008 Experimental data were obtained outside the VCO clipping regions and for values of the delay for which the in-phase synchronized state is stable (σ < 0). Fig 11 shows the experimental and theoretical results for the perturbation response rate σ. Again, the experimental results are in quantitative agreement with the results from the phase model.

Ring of three mutually coupled DPLLs
We now study a ring of three mutually coupled DPLLs (no. 1, 2, and 3 in Table 1). In this experimental setup, we observe in-phase, 1-twist, and equivalent 2-twist synchronized states. Fig 12 shows the measured collective frequency O of the in-phase (0-twist) and the 1-twist state for different values of the transmission delay, together with the results from the phase model, Eq (20). The behavior of these solutions is in quantitative agreement with the results of the phase model.

Square lattice of 3 × 3 mutually coupled DPLLs
We now study a square lattice of 3 × 3 mutually coupled DPLLs with periodic and open boundary conditions (no. 1-9 in Table 1). Depending on the type of boundary condition, we observe in-phase, m-twist, or checkerboard synchronized states. For a nearest-neighbor coupled system in two dimensions, m-twists can occur in each dimension individually, corresponding to constant phase offsets θ m 1 and θ m 2 in x and y-direction, respectively. For the special case of θ m 1 = θ m 2 , the results for the collective frequency O are identical to those for a one-dimensional ring discussed in Sec. 3. Similarly, the checkerboard synchronized state in two dimensions has the same collective frequency as the checkerboard state in one dimension.

Conclusion and discussion
In this paper, we have shown in theory and in experiments that global synchronization can be achieved by mutually coupling DPLLs with transmission delays. We used a phase oscillator description for networks of digital PLLs, taking into account the effects of filtering and transmission delays, to analyze frequency-synchronized states. The collective frequencies of inphase, m-twist, and checkerboard synchronized states were obtained analytically. We analyzed the stability of those states and showed, how the DPLL's filters and transmission delays in the network affect synchronization, in particular the collective frequency and the time scales of synchronization. Specifically, the collective frequency of a synchronized state in general differs from the intrinsic frequency of the DPLLs. In the presence of transmission delays the collective frequency of a synchronized state deviates from the intrinsic frequencies of the DPLLs. We found that it is essential to take into account the dynamic properties of filtering as it has profound effects on the stability properties. In particular, filtering introduces new stability transitions compared to delay-coupled phase oscillators without filters, which were discussed in [38]. Filtering leads to an inert system behavior. This can most easily be seen for the case of filters of first order, a = 1, in which the dynamic equations can be expressed by a second-order differential equation comprising an inertial term, see C. For spatially extended systems of coupled phase oscillators with inertia, it has been shown that perturbations can initiate transient waves that travel through the system and that might affect the synchronization behavior [31]. Our approach of mutually coupled distributed DPLLs enables frequency-synchronized states for many oscillators and thus has the potential to scale up to large systems. It is particularly relevant for, e.g., massive MIMO systems due to considerable delays induced by typical spacing between antennas. However, it can also be applied to other systems where a large number of clocks has to be synchronized.
We designed experiments to test the predictions of our phase model in real systems of coupled DPLLs. Our measurements of the collective frequency and the synchronization relaxation times show that our theory can quantitatively describe the behavior of mutually delay-coupled DPLL systems. All dynamic states discussed here in theory have been shown to exist in real networks of coupled DPLL systems in our experiments. Centralized state-of-the-art systems require a master clock together with complicated clock trees to ensure that the entrained DPLLs receive the clock signal at the same time. In contrast, the system proposed here relies on a regular lattice geometry and is therefore readily implemented and trivially scalable. Moreover, wiring distances between components are kept to a minimum. Through systematically taking into account the effects of transmission delays that enable synchronization, these delays become a design parameter rather than a limitation of the system.
The theory presented here is a practical tool to design systems of coupled DPLLs for technical applications and to predict their behaviors. We showed that with XOR phase detectors, the stability of such synchronized states depends strongly on the transmission delays in the coupling and the properties of the loop filter. Other theoretical studies for networks of mutuallycoupled PLLs with flip-flop phase detectors report that the effects of time-delays can be  neglected [44]. We have considered an idealized system of delay-coupled DPLLs without phase noise. Phase noise is present in any real system of coupled DPLLs and can lead to fluctuations such as cycle slips. While our experiments, in which phase noise is inevitable, are nevertheless well described by our theory, a detailed study of the effects of phase noise is currently conducted in our group. The problem of efficiently booting such systems into synchrony and the role of phase noise are interesting topics for future research.

A Approximation of the control signal
Here we determine the frequency contributions of the PD's output signal x PD k ðtÞ and discuss the impact of the low pass filter, i.e. the damping, on the these components. The Fourier representation of the square wave function P and the triangle function Δ are given by Eqs (2) and (7). Hence, evaluating Eq (3) yields,  Self-organized synchronization of delay-coupled DPLLs where a i = 2i + 1 and ϕ l,τ (t) = ϕ l (t − τ). The different sums in the second identity represent different frequency components. From Eq (5), it can be seen that _ 0 k À _ 0 l ¼ OðK VCO Þ. Hence, the first sum in the second identity contains low frequency components. The lowest frequency components in the second sum are given by the terms with i + j = 1. From Eq (5), it can be seen that a i _ 0 k À a 1À i _ 0 l ¼ Oð2o 0 Þ for i 2 {0, 1} and hence, the second sum contains only high frequency components. The lowest frequency component in the third sum is given by the term with i = j = 0. From Eq (5), it can be seen that _ 0 k þ _ 0 l ¼ Oð2o 0 Þ and hence, the third sum contains only high frequency components. Thus, using Eq (7), we can write x PD k as x PD k ðtÞ ¼ where R HF denotes the high frequency components of the signal. In Eq (6), we approximate R 1 0 du pðuÞ R HF ðt À uÞ ' 0, corresponding to ideal suppression of the high frequency components by the loop filter [32].

B Collective frequency of the m-twist synchronized state
We state the explicit solutions to Eq (20) for the collective frequency O of the m-twist synchronized states. We use the piecewise linear behavior of the rhs of Eq (20) to obtain these solutions as intersection points between the lhs and rhs of Eq (20), analogously to the solutions for the in-phase synchronized state given in Section 2.1.
We define the constants μ and ν by Potential intersection points between the lhs and the rhs of Eq (20)  O À Ã is a solution if and only if the condition P À p; n½1 À 2my m þ 2mn À 1 þ n 2 holds. O AE j is a solution if j falls within a range j 2 N between A ν and B ν − 1/2 for O þ j and between A ν and B ν + 1/2 for O À j , where A ν = (ω − νK)τ/2π and B ν = (ω + νK)τ/2π and, in addition, the condition Pðy m À p; 2mp À y m ; O AE j tÞ ð38Þ holds.

C First order signal filtering in a PLL and inertia
Models of coupled oscillators with inertia are an active and emerging research field [45][46][47][48].
Here we show that for a LF of first order, the filter integral in the phase model Eq (8) can be rewritten and yields a second order delay-differential equation including inertia. The phase model for first order filters with a = 1, given in Eq (8), and finite phase history starting at t = 0 can be expressed in Laplace space as where0 k ðlÞ denotes the Laplace transform of the phase andx PD k ðlÞ the Laplace transform of the phase detector signal. Transformation back into time domain yields b € 0 k ðtÞ þ _ 0 k ðtÞ ¼ o þ ½o À _ 0 k ð0Þbdðt À 0Þ þ Kx PD k ðtÞ ; ð40Þ where x PD k ðtÞ ¼ n À 1 k P l d kl h½0 l ðt À tÞ À 0 k ðtÞ. The term ½o À _ 0 k ð0Þbdðt À 0Þ denotes a phase kick at t = 0 with b, the characteristic integration time of the filter, which shifts the phases of all oscillators determined by the initial conditions. Comparison with Eq (8) reveals that the signal filtering with a filter of first order can be rewritten as an inertia term of the phase variable.