Onset of synchronization in networks of second-order Kuramoto oscillators with delayed coupling: Exact results and application to phase-locked loops

We consider the inertial Kuramoto model of $N$ globally coupled oscillators characterized by both their phase and angular velocity, in which there is a time delay in the interaction between the oscillators. Besides the academic interest, we show that the model can be related to a network of phase-locked loops widely used in electronic circuits for generating a stable frequency at multiples of an input frequency. We study the model for a generic choice of the natural frequency distribution of the oscillators, to elucidate how a synchronized phase bifurcates from an incoherent phase as the coupling constant between the oscillators is tuned. We show that in contrast to the case with no delay, here the system in the stationary state may exhibit either a subcritical or a supercritical bifurcation between a synchronized and an incoherent phase, which is dictated by the value of the delay present in the interaction and the precise value of inertia of the oscillators. Our theoretical analysis, performed in the limit $N \to \infty$, is based on an unstable manifold expansion in the vicinity of the bifurcation, which we apply to the kinetic equation satisfied by the single-oscillator distribution function. We check our results by performing direct numerical integration of the dynamics for large $N$, and highlight the subtleties arising from having a finite number of oscillators.


I. INTRODUCTION
A. The model The Kuramoto model with inertia is representative of complex many-body dynamics involving a set of rotors characterized by their phases and angular velocities that are coupled all-to-all through the sine of their phase differences. Specifically, the dynamics for a system of N rotors is given by a set of 2N coupled first-order differential equations of the form [1][2][3] where the dot denotes derivative with respect to time, θ i ∈ [0, 2π) and v i are the phase and the angular velocity of the i-th rotor, respectively, whose moment of inertia is m > 0. Here, γ > 0 is the damping constant, K > 0 is the coupling constant, while ω i ∈ [−∞, ∞] is the natural frequency of the i-th rotor. The frequencies {ω i } 1≤i≤N constitute a set of independent and quenched disordered random variables distributed according to a given distribution g(ω), normalized as of the Kuramoto model are N coupled first-order differential equations of the form The mean-field nature of either the dynamics (2) or the dynamics (1) becomes evident on defining the so-called synchronization order parameter R(t) and the global phase Φ(t), as [4] with 0 < R < 1 characterizing a synchronized phase, and R = 0 an incoherent phase. In terms of R(t), the dynamics (1) may be rewritten asθ mv i (t) = −γv i (t) + γω i + KR(t) sin(Φ(t) − θ i (t)), which shows that the evolution of the dynamical variables at time t is governed by the value of the mean-field R(t)e iΦ(t) set up collectively at time t by all the oscillators.
Both the models (1) and (2) have been extensively studied in the past and a host of results have been derived with regard to the parameter regimes allowing for the emergence of a synchronized stationary state (see Ref. [10] for a recent overview). For example, consider a g(ω) that is unimodal, namely, it is symmetric about its mean ω 0 , and decreases monotonically and continuously to zero with increasing |ω − ω 0 |. In this case, it is known that in the stationary state of the dynamics (2), the system for a given choice of g(ω) may exist in either a synchronized or an incoherent phase depending on whether the coupling K is respectively above or below a critical value K c = 2/(πg(ω 0 )); on tuning K across K c from high to low values, one observes a continuous phase transition in R st , the stationary value of R(t). Namely, R st decreases continuously from the value of unity, achieved as K → ∞, to the value zero at K = K c and remains zero at smaller K values. One may interpret the transition as the case of a supercritical bifurcation, in which on tuning K, a synchronized phase bifurcates from the incoherent phase at K = K c . In particular, a small change of K across K c results in only a small change in the value of R st ∝ √ K − K c close to and above K c [4,16]. For the same choice of a unimodal g(ω), the inertial dynamics (1) on the other hand show a discontinuous phase transition between synchronized and incoherent phase, where R st exhibits an abrupt and big change from zero to a non-zero value on changing K by a small amount across the phase transition point [17,18]. Here, the bifurcation of the synchronized from the incoherent phase is said to be subcritical and leads to hysteresis [19]. Thus, presence of inertia is rather drastic in that it changes completely the nature of the bifurcation and hence of the underlying stationary state.
In this work, we study for the first time the effect of a delay in the interaction between the oscillators within the framework of dynamics (1). The dynamical equations of this modified model are given bẏ thereby modeling the time-evolution which is governed by the value of the mean field at an earlier instant t − τ , where τ > 0 is the time delay in the interaction between the oscillators. Here, α ∈ (−π/2, π/2) is the so-called phase frustration parameter, an additional dynamical parameter that is known to affect significantly the behavior of the Kuramoto model [20]. In the overdamped limit, the model (5) reduces to which in presence of additional Gaussian, white noise has been addressed in Ref. [21]. Note that for the dynamics in (6), the parameter γ may be scaled out by a redefinition of time, so that the relevant dynamical parameters are K, τ and α. In a recent work [22], two of the present authors have investigated the dynamics (6), deriving for generic g(ω) and as a function of the delay exact results for the stability boundary K c (τ ) between the incoherent and the synchronized phase and the nature in which the latter bifurcates from the former at the phase transition point (the effects of changing τ at a fixed α are the same as those from changing τ at a fixed α keeping α − ω 0 τ constant). Our results imply that for a given choice of g(ω), the nature of transition (continuous versus discontinuous) between the synchronized and incoherent phases depends explicitly on the value of τ .
In view of the aforementioned developments, it is evidently of interest to investigate the effects of inertia on the time delayed model and thus embark on a detailed analysis of the dynamics (5). Since even without delay, inertia is known to have nontrivial and interesting consequences as mentioned above, we may already anticipate that an interplay of the influence of delay and inertia may result in an even richer stationary state for the dynamics (5) vis-à-vis for dynamics (6). Remarkably, the dynamics in Eq. (5), far from being just a model of academic interest, emerge naturally in the context of coupled phase-locked loops, as we now demonstrate.
B. Relation to a network of phase-locked loops A phase-locked loop (PLL) is an electronic component designed to generate an output signal that has a constant phase relation (and is thus locked) to the phase of its input reference. Fig. 1 shows a schematic PLL architecture consisting of a phase detector (PD), a loop filter (LF), and a voltage-controlled oscillator (VCO) acting as a variablefrequency oscillator, all connected in a feedback loop. The phase-detector output x PD (t) represents the phase relations of the periodic output signal x out (t), generated by the VCO, with the phase of the periodic input-signal x in (t). The loop-filtered phase-detector output yields the control signal x C (t) that controls the instantaneous frequency of the VCO so that its corresponding output approaches the phase and frequency of the input signal. The latter property enables a PLL to track an input frequency, or, to generate a frequency that is a multiple of the input frequency. PLL's find wide use in electronic applications as an effective device to, e.g., recover a signal from a noisy communication channel, generate a stable frequency at multiples of an input frequency, and to distribute a quartz reference clock signal via a clocktree architecture.
Let us now consider the setup of N ≥ 2 mutually delay-coupled PLL's occupying the nodes of a network, in which the input signals for a given PLL are constituted by the delayed output received from other PLL's [23][24][25]. The delay could be due to transmission signaling-times, and is accounted for in the following by a discrete delay-time τ . We consider the LF to ideally damp the high-frequency components of the PD signal. Consider the output signal of the i-th PLL, i = 1, 2, . . . , N , in the network where θ i (t) denotes the phase, and s(θ i (t)) is a 2π-periodic function with amplitude one. Depending on the type of PLL, i.e., analog or digital, the output signal may be sinusoidal or a rectangular function, respectively. The VCO is operated such that its output frequencyθ i (t) depends linearly on the control signal x C i (t): where ω VCO i,0 denotes the natural frequency, and K VCO i the VCO input sensitivity. The control signal is the output of the loop-filter: where x PD i (t) denotes the phase-detector signal, and p(u) is the impulse response of the filter. Considering first-order loop-filters, i.e., p(u; 1, b) being the Γ-distribution with shape parameter a = 1 and scale parameter b, the above integral equation can be rewritten by using Laplace transforms [24,26], yieldinġ where ω c denotes the cut-off frequency of the first-order low-pass filter. The initial state of the filter is given by The phase-detector signal depends on the type of PLL: where C is a PLL type specific offset (C = 1/2 for XOR PD's, while C = 0 for multiplier PD's), c ij = {0, 1} are the components of the adjacency matrix, with the value 1 (respectively, 0) denoting whether PLL units i and j are coupled (respectively, uncoupled), n(i) ≡ j c ij the total number of units coupled with unit i, h(x) is a 2π-periodic coupling function, and we assumed the high-frequency components to be filtered ideally by the LF [23]. Equations (8)- (11) combined together yield a second-order phase model with delayed-coupling: where we have defined ω i ≡ ω VCO i,0 The 2π-periodic coupling function h depends on the type of the PD and the corresponding input signals. Here we consider the case of a cosine coupling function, h(x) = cos(x), for analog PLL's and multiplier phase-detectors and a triangular coupling function, and h(x) = ∆(x) for digital PLL's with XOR phase-detectors. In the latter case the coupling-function can be approximated as h(x) = −8/π 2 cos(x). The case of a d-flip flop [36] phase detector for digital PLLs, which has a linear coupling-function, will not be considered in this work. Given these cases, we will use a sinusoidal coupling-function with a phase frustration parameter α ∈ [−π/2, π/2], that is, with h(x) = sin(x − α), which represents both of the cases mentioned above. We will also specialize to the case when every PLL unit is coupled to every other, implying c ij = 1 ∀ i, j = 1, 2, . . . , N and n(i) = N . Comparing Eqs. (12) and (5) leads to the correspondence m = ω −1 c , γ = 1, as well as K = K i , α = π/2 for the analog PLL case, and α = −π/2, K = 8 K i /π 2 for the digital PLL approximation.
Before moving on to an analysis of the dynamics (5), it is pertinent that we give here a summary of our results obtained in this paper and the techniques employed in achieving them. We here obtain exact analytical relations for the critical point K c (τ ) beyond which the incoherent phase of the dynamics (5) becomes unstable, and furthermore, the nature in which the synchronized phase bifurcates from the incoherent phase as K is increased beyond K c (τ ). An illustration of our results for a unimodal Lorentzian distribution is shown in Fig. 2 for two representative values of the inertia, which displays both K c and the real part of a quantity c 3 whose sign determines the nature of the bifurcation of the order parameter R, Eq. (3): a positive (respectively, a negative) sign implies a subcritical bifurcation and hence, a discontinuous transition (respective, a supercritical bifurcation and hence a continuous transition). As may be seen from the Fig. 2, K c and the real part of c 3 both have an essential dependence on τ and m, while our analysis (see Eq. (21)) suggests that the effects of changing τ at a fixed α are the same as those from changing τ at a fixed α keeping α − ω 0 τ constant.
We now summarize our method of analysis in obtaining the aforementioned results. We start off with considering the dynamics (5) in the limit N → ∞, when it may be effectively characterized by a single-oscillator probability density F (θ, v, ω, t), which gives at time t and for each ω the fraction of oscillators with phase θ and angular velocity v. The time evolution of F (θ, v, ω, t) follows a kinetic equation, of which the incoherent state f 0 (θ, v, ω) (associated with R st = 0) represents a stationary solution. We rewrite the kinetic equation in the form of a delay differential equation (DDE) [27,28] for perturbations f t (ϕ) ≡ F (θ, v, ω, t+φ); −τ ≤ φ < 0 around f 0 (θ, v, ω). The DDE involves a linear evolution operator D and a nonlinear one, F . We obtain the eigenvalues and the eigenvectors of D and of the corresponding adjoint operator D † . As is well known [19], the knowledge of the eigenvalues allows to locate the critical value K c of the coupling K above which the incoherent state f 0 (θ, v, ω) becomes linearly unstable. We then build for K > K c the unstable manifold expansion of the perturbation f t (φ) along the two complex conjugated eigenvectors associated with the instability. Using a convenient Fourier expansion of the relevant quantities and working at K slightly greater than K c , we thus obtain the amplitude dynamics describing the evolution of perturbations f t (ϕ) in the regime of weak linear instability, K → K + c . The nature of the amplitude dynamics at once dictates the nature of bifurcation occurring as soon as K is increased beyond K c : The amplitude dynamics has a leading linear term and a nonlinear (cubic) term, and as is well known from the theory of bifurcation [19], the sign of the real part of this cubic term (denoted c 3 in Fig. 2) dictates the nature of the bifurcation, with positive and negative signs leading respectively to subcritical and supercritical bifurcation.
The paper is organized as follows. Section II forms the core of the paper, in which we derive our main results, Eqs. (18) and (21), with most technical details of the computation relegated to the four appendices. We illustrate our analytical results with the representative example of a unimodal Lorentzian distribution. In Section III, we make a detailed comparison of our analytical results obtained in the limit N → ∞ with numerical results for finite N obtained by performing numerical integration of the equations of motion. Here, in particular, we discuss the subtleties involved in making such a comparison whose origin may be traced to finite-size effects prevalent for finite N . The paper ends with conclusions.

II. EXACT ANALYSIS IN THE LIMIT N → ∞
We now turn to a derivation of our results for the system (5). To this end, consider the system in the limit N → ∞, when the dynamics may be effectively characterized in terms of the single-oscillator probability density F (θ, v, ω, t) defined above. This density is 2π-periodic in θ, and obeys the normalization The time evolution of F (t) ≡ F (θ, v, ω, t) may be derived by following the procedure given in Ref. [7]. One obtains the evolution equation (14) where we have defined as functionals of F the quantity In particular, R −1 coincides with the N → ∞ limit of the Kuramoto order parameter (3). From Eq. (14), one may check that the incoherent state solves the equation in the stationary state and thus represents an incoherent stationary state. To examine how in the stationary state the incoherent stable becomes unstable as K is tuned above a critical value K c , we employ an unstable manifold expansion of perturbations about the incoherent state in the vicinity of the bifurcation. To perform the analysis, we write F = f 0 + f , with f being the perturbation, and rewrite Eq. (14) as a sum of a part D linear in f and a nonlinear part F : Here, we have defined f t (φ) ≡ f (t + φ), while the expressions of the two operators D and F are given in Appendix VI A. Starting with Eq. (17), the unstable manifold expansion involves a linear and a weakly nonlinear analysis, and requires combining two formalisms: i) the one developed in Ref. [18] for the case of the Kuramoto model with inertia but with no delay, ii) the delay formalism [27][28][29], as done in Ref. [22]. The steps are detailed in the Appendices. We summarize here the main results. The linear stability analysis yields the dispersion relation (see Appendix VI C) which has its roots giving the eigenvalues associated with the linear operator D. In particular, for K ≥ K c , the stationary state f 0 becomes unstable, with associated unstable eigenvalues λ satisfying Re(λ) ≥ 0. Note that for K < K c , the incoherent state is neutrally stable, i.e., there is no discrete eigenvalue but only a continuous spectrum; in this case, perturbations f 0 are damped in time via a mechanism similar to the Landau damping [30]. The weakly nonlinear analysis describes the type of bifurcation as K → K + c and hence as Re(λ) → 0 + . The analysis involves decomposing the perturbation into a contribution along the unstable eigenvectors P (θ, v, ω, φ), P * (θ, v, ω, φ) associated with the unstable eigenvalues λ, λ * and a contribution S t (θ, v, ω, φ) in the perpendicular direction, as where c.c. stands for complex conjugation, A(t) = (Q, f t ) τ is the amplitude of the unstable mode, and (Q, S t ) τ = 0.
Here, we have introduced the eigenvector Q of the adjoint operator D † and the scalar product (·, ·) τ , see Appendix VI B. The unstable manifold approach consists in expanding the perpendicular component S t in terms of the small amplitude A, S t (θ, v, ω, φ) = H[A, A * ](θ, v, ω, φ) and computing H perturbatively. For small A (that is, in the close vicinity of K c ), it can be shown that R(t) = A * (t) + O(|A| 2 A * (t)), so that studying the bifurcation of A is equivalent to that of the order parameter R. Following the nonlinear study based on ideas developed in Refs. [18,22] and detailed in Appendix VI D, we find the following reduced equation for the order parameter: where the unstable eigenvalue λ is decomposed into its real and imaginary parts: λ = λ r + iλ i . A few remarks are in order: a) The coefficient c 3 diverges as λ r → 0, which is the regime where the reduction is valid. This singular behavior is typical of this type of systems [16,18,31], and stems from the existence of the continuous eigenspectrum that cannot be described by the finite dimensional equation (20). b) However, we still expect the behavior of c 3 to determine the type of bifurcation. For Re(c 3 ) > 0, we expect a subcritical (discontinuous) bifurcation, while for Re(c 3 ) < 0, we expect a supercritical bifurcation. In the latter case, the scaling of the stationary amplitude is A st ∝ λ r , which differs from the usual Kuramoto model where it goes as √ λ r . c) In the case with no inertia, that is, with m = 0, we expect the coefficient c 3 to be quantitatively relevant in giving the exact amplitude A st of the stationary branch close to the bifurcation; here, because of the singularity, only the sign and scaling of c 3 (λ) can be used to get qualitative information.
A. Application to a Lorentzian distribution To assess the effects of inertia in a delay system, we consider a Lorentzian distribution of the natural frequencies: g L (ω) = σ/[π((ω − ω 0 ) 2 + σ 2 )]. The dispersion relation (18) at criticality gives K = K c and λ = 0 + + iλ i,c : Solving this system gives us K c (τ ) and λ i,c (τ ). Then the sign of the cubic coefficient Re(c 3 )(τ ) is given by where we used Eq. (21) with g = g L . We plot in Fig. 2 the quantities K c and Re(c 3 ) as a function of the delay for two different inertia values m = 0.1 and m = 3 s/rad for a Lorentzian distribution. For τ → 0, we recover the no-delay results showing a positive Re(c 3 ) and hence a subcritical bifurcation [7]. Moreover, as in the case with no inertia [22], m = 0, the delay induces "oscillations" in the sign of Re(c 3 ). We observe that different nonzero values of m do not change much the behavior of the bifurcation.

A. Method
In the preceding section, we provided an analytic characterization of the stability properties of the incoherent state in the Kuramoto model with delayed coupling and inertia in the limit N → ∞. Here we present results from numerical integration of the dynamics (5) with Lorentzian-distributed natural frequencies with location parameter ω 0 = 3 radHz and scale parameter σ = 0.1 radHz. For all-to-all coupling, Eq. (5) can be rewritten in terms of the Kuramoto order parameter: Using R x (t − τ ) = 1/N j cos(θ j (t − τ )) and R y (t − τ ) = 1/N j sin(θ j (t − τ )), we rewrite the equations of motion asθ As mentioned earlier, the effect of α for a fixed τ can be studied by considering α − ω 0 τ instead of τ in the dynamics, hence here we have set α to zero without loss of generality. Hence, for γ = 1, we have the set of equationṡ which we integrate numerically using an Euler iteration-method, given the initial phases θ i (0) independently and identically distributed in [0, 2π), and initial state of the filters x C i (0). R x (t hist ) and R y (t hist ) with t hist ∈ [−τ, 0] are given by the history of the network of oscillators, which we obtain by evolving each oscillator independently according to its own natural frequency. The code is written in python and compiled with cython for fast execution and can be found on the Gitlab repository here [32].
Our objective behind performing the numerics is to verify our theoretical results obtained in the thermodynamic limit for the critical coupling constant K c above which the incoherent state becomes unstable. Furthermore, we want to confirm the type of bifurcation as predicted by our theoretical results that can be observed as the coupling constant K is tuned across K c . To this end, we integrate numerically the dynamics of large networks of all-to-all delay-coupled PLL clocks in the vicinity of the theoretically-predicted K c , see Fig. 2.
We proceed with our numerical work as follows. For a given set of parameters (N, τ, m, γ = 1, α = 0), a vector of discrete coupling constants K n = {K on , K on + ∆K, . . . , K end − ∆K, K end }; ∆K > 0 and a set of Lorentziandistributed natural frequencies {ω i }, each oscillator is evolved independently with its own natural frequency for a time τ to obtain the dynamical history for N PLLs in the interval [−τ, 0]. In the next step, we turn on the coupling between the oscillators at an initial coupling constant K on that is close to but smaller than the critical coupling constant predicted by our theoretical results. Subsequently, the delay-coupled system of all-to-all coupled PLLs is evolved with the coupling constant kept fixed at K on for time T num that is long compared to the mean period of the independent oscillators, in order to ensure that the system settles into a stationary state at the fixed value of the coupling. Then, using the phases of the final interval [T num − τ, T sim ] as the history, we evolve the system of coupled oscillators for the next larger value in K n for time T num , and so on, until the value K end is reached. In the final part of this procedure, we follow the exact reverse protocol, namely, repeating the above steps while decreasing the value of the coupling from K end until the value K on < K c is reached. In numerics, we track the value of the Kuramoto order parameter in time, and save for each value of the coupling in the set K n the final value of the order parameter obtained at the end of run for time T num as well as its average and variance computed over a time equal to 50 times the time period corresponding to ω 0 .  In our simulations we find the bifurcations that were predicted by the theoretical results. As the coupling constant K increases, it crosses a critical coupling constant K num c as found in our finite-size simulation, and we observe a subcritical (discontinuous) transition and hysteresis for Re(c 3 ) > 0. For Re(c 3 ) < 0 on the other hand, we find a supercritical (continuous) transition with a linearly growing order parameter and no hysteresis as K grows larger than K num c . We observe that for the case of m = 0.1 s/rad, the hysteresis seems to be weaker than in the case with m = 3.0 s/rad. The results are in good agreement with our theoretical predictions for K c and the type of bifurcation.

C. Discussion of numerical results
Numerical validation of our analytical results obtained in the thermodynamic limit is not trivial and comes with a few difficulties that we will discuss here. Since we do not know the system size N crit of oscillators below which strong finite-size effects will come into play and the number N thermo above which the behavior coincides with that in the thermodynamic limit, we decided while doing our numerical checks to go for as large a system size as is practicable. This however becomes very resource intensive, since for systems with time delays, a memory of the states of all oscillators for a time period [t − τ, t] has to be stored in order to perform dynamical evolution. For the mean-field coupling case, we have the advantage that it is sufficient to keep only the history of the order parameter variables R x (t − τ ) and R y (t − τ ).
It is known from the literature [33] that for the Kuramoto model in absence of delay and inertia, when the number of oscillators is finite, we may expect to find K num c ≤ K theory c , depending on the number N of oscillators considered. This difference is even stronger when considering the inertial model and subcritical bifurcation, where the convergence [17] is very slow with N (compared with K theory c (N ) − K num c (N ) ∼ N −0.4 without inertia [33]). Note that the prediction K theory c > K num c was obtained for systems without delay, thus observing K theory c ≈ K num c in Fig. 3 and Fig. 4 at τ = 2 s does not contradict the results in [17,33] and raises an interesting issue for future investigation as to how the difference between K num c and K theory c scales with N . However, in cases with supercritical bifurcation, the numerically observed value of the critical coupling seems more different from the theoretical value than in cases with subcritical bifurcation (which could be another indicator of the type of transition). We thus have no prior knowledge of the exact value K num c (N ) at which to expect the transition. Furthermore, the multistability present in the system due to the delayed interaction results in a number of step-like transitions that follow once the incoherent state becomes unstable. As a consequence, validation of a linear and continuous transition for K > K c becomes hard to resolve when increasing K n further than the regime of the continuous transition. For smaller values of the inertia, the discontinuous transitions and hysteresis regimes also become much smaller, and it becomes difficult to resolve them even with, e.g., N = 10 5 oscillators. Furthermore, there are no a priori conditions to guide us while choosing the values of ∆K and T sim . The tests required in estimating the right choices for the dynamical parameters are computationally expensive and require long computation times.

IV. CONCLUSIONS
In this work, we have studied the effect of time delay in the interaction between oscillators within the framework of the inertial Kuramoto model of globally coupled oscillators. For a generic choice of the natural frequency distribution of the oscillators, we obtain exact analytical results that imply that in contrast to the case with no delay, the system in the stationary state may exhibit either a subcritical or a supercritical bifurcation between a synchronized and an incoherent phase. The precise nature of bifurcation has an essential dependence on the amount of delay present in the interaction as also on the value of inertia of the oscillators. Our theoretical analysis, performed in the limit of an infinite number of oscillators, is carried out by employing an unstable manifold expansion in the vicinity of the bifurcation, which we apply to the kinetic equation satisfied by the single-oscillator distribution function, Eq. (14). The one-dimensional reduction, Eq. (20), of the dynamics for the order parameter is plagued by singularities that are reminiscent of an infinite dimensional bifurcation and, thus, gives at best qualitative information on the bifurcation nature, a fact that our numerical results fully support. We notice, however, that the unstable manifold method is very robust in the context of kinetic equations with continuous spectrum, since it is, to the best of the authors' knowledge, the only one giving analytic predictions for the Kuramoto model both with inertia m = 0 and delay τ = 0, while other methods, like the Ott-Antonsen ansatz [34], work only for m = 0, while self-consistent methods e.g. [1] have been applied only for the case with no delay, τ = 0. Direct numerical integration of the dynamics allows to highlight the subtleties one is confronted with when checking the analytical results against those obtained numerically for a finite number of oscillators. For systems of delay-coupled PLLs with heterogeneous natural frequencies, our results allow to predict the minimal coupling sensitivity of the voltage-controlled oscillators necessary to enable the network to become synchronized. Moreover, such PLL networks generally seem to exit the incoherent state at smaller coupling sensitivity if the transition happens through a subcritical bifurcation, Re(c 3 ) > 0, and close to integer multiples of the mean natural period of the oscillators, where we find the local minima of K c , see Fig. 2. It may be noted that with increasing transmission delay, the onset of synchronization can generally be achieved at smaller values of K c . Hence, larger values of the transmission delay seem to decrease the stability of the incoherent state. The time evolution of a function F (t) according to a nonlinear operator with delay, M [F (t)], can be rewritten in term of a delay variable ϕ such that the time-evolution operator is given by with F t (φ) ≡ F (t + φ). Expanding around a stationary state f 0 , as F = f 0 + f , with f being a perturbation, we define the linear and nonlinear operators D and F , according to We decompose the linear operator into two parts, L = L + R, namely, a part L that does not contain any delay term and a part R that has all the delay terms. Rewriting Eq. (14) according to the above formalism yields where we use the shorthand ∂ v ≡ ∂/∂v for derivatives and use in all appendices the transformation m → m/γ and K → K/γ. Moreover, to simplify the Appendices, we work in the rotating frame θ j (t) → θ j (t)−ω 0 t, ω j → ω j −ω 0 ∀ j, so that the distribution g(ω) is now centerer in 0.

B. Dual space
In the functional space of delayed functions, there is no L 2 canonical inner product. However, Ref. [27] defines a bilinear form acting as the inner product on this space. In our problem with a discrete delay, the scalar product is where (q(0), p(0)) denotes the usual scalar product on L 2 (T × R × R) (phase, angular velocity and natural frequency) and the integral term contains the delay contribution. The adjoint of the linear operator D, obtained by using the equality (q(ϕ), Dp(ϕ)) τ = (D † q(ϕ), p(ϕ)) τ , is defined in the dual space, and is given by We also decompose L † = L † + R † , with

D. Unstable manifold Expansion
To perform the weakly nonlinear analysis, we will expand the perturbation f t (φ) along the unstable eigenvector and the unstable manifold. The decomposition on the unstable manifold reads with A(t) = (Q, f t ) τ , (Q, P * ) = 0 and (Q, H) = 0. We assume that H is at least of order (A, A * ) 2 . The amplitude A is directly related to the Kuramoto order parameter R as R = A * + O(|A| 2 A * ). Let us define the following Fourier expansions needed for further analysis: where the dependance on A of the Fourier coefficients of H is imposed by rotational symmetry [31]. To proceed with the analysis, we will need to expand the coefficients w k in powers of |A| 2 , w k = ∞ j=0 |A| 2j w k,j . To be consistent with the assumption of the unstable manifold being at least of order (A, A * ) 2 , we need to have w ±1,0 = 0.
The Fourier coefficients of the nonlinear operator (34) are Note that contrary to the case with no inertia, m = 0, where L 0 = N 0 = 0 so that (f t ) 0 = constant = 0, it is not so in the present case so that (f t ) 0 = 0 and w 0 = 0. This difference will have major consequences for the reduction, like a 1/λ divergence in the c 3 coefficient. For ϕ = 0, we find w 0 (ϕ) = h 0,0 e 2λrϕ + O(|A| 2 ), w 2,0 (ϕ) = h 2,0 e 2λϕ + O(|A| 2 ), and with the boundary equation ϕ = 0: Solving these equation will give us h 0,0 and h 2,0 needed in the following. Projection of the dynamics along the unstable mode using (Q, (31)) τ yields the equation for the amplitude A(t) to bė where we used Eq. (51) for k = 1 keeping only the leading order. To determine the nature of the bifurcation, we must compute explicitly the coefficient c 3 . To do that, we must first compute the Fourier component of the unstable manifold.

Putting everything together
One can ascertain that the only diverging term will come from ψ (2) * W * 1 dω; thus, the leading term is We conclude that the leading behavior of c 3 for m > 0 is given by In particular, the sign of s(τ ) ≡ Re e i(α−(ω 0 +λ i )τ )