Delay master stability of inertial oscillator networks

Time lags occur in a vast range of real-world dynamical systems due to finite reaction times or propagation speeds. Here we derive an analytical approach to determine the asymptotic stability of synchronous states in networks of coupled inertial oscillators with constant delay. Building on the master stability formalism, our technique provides necessary and sufficient delay master stability conditions. We apply it to two classes of potential future power grids, where processing delays in control dynamics will likely pose a challenge as renewable energies proliferate. Distinguishing between phase and frequency delay, our method offers an insight into how bifurcation points depend on the network topology of these system designs.

Introduction. The study of nonlinear dynamics on complex networks has received full interdisciplinary attention in past years due to its potential for modeling the complexity of real-world dynamical systems. An intrinsic feature of such systems is that their time evolution generally depends on past states. Time delays, caused by finite propagation speeds or processing times, induce retarded reactions of variables to changes in the system. For example, delays occur in lasers because of the finite speed of light [1,2]; population dynamics depend on maturation and gestation times [3], and the exchange of information between neurons requires time for both signal transmission as well as processing [4].
Mathematically, continuous delay problems are described by delay differential equations (DDEs) [5,6]. From their analysis it is known that delays can substantially alter a system's asymptotic behavior [7]. However, asymptotic stability analysis of DDEs is challenging because the corresponding spectrum contains an infinite number of complex roots. In fact, exact conditions for stability pose an open problem in research, especially regarding networks. Most previous studies have been limited to numerical investigations of characteristic equations or restricted to simple network topologies, often yielding only sufficient stability criteria.
Recent work has led to a thorough analytical understanding of the spectrum in the limit of large delay, with applications in e.g. optoelectronics [8][9][10][11]. In many cases, however, time lags may match the system's dynamical timescales and may play a critical role for stability. Particularly in systems of coupled oscillators like the paradigmatic Kuramoto model [12], delays often become comparable to the oscillation period. There, the asymptotic stability of a synchronous regime is a central property with crucial implications for applications.
Pecora and Carroll have developed a powerful method known as the master stability formalism to determine the * Email: hellmann@pik-potsdam.de stability of synchronization for identical oscillators without delay [13]. The main idea is to project the state vector into the eigenspace of the coupling matrix, yielding a block diagonal form that defines the associated master stability function (MSF). This way, dynamical parameters of the system are separated from topological information about the network.
Several studies have calculated MSFs for specific models with time-delayed couplings [14][15][16][17]. Here, for the first time, we generalize the formalism to DDE inertial oscillator models containing an arbitrary constant discrete delay τ > 0 that may appear in the local dynamics as well as in a diffusive coupling term. While the master stability formalism requires complete synchronization of oscillators, we merely assume phase synchronization where oscillators may have constant phase differences [18]. Our analytic approach leads to a decomposition into second-order DDEs in terms of the eigenvalues of the graph Laplacian matrix. Based on results from Bhatt and Hsu [19], we derive necessary and sufficient conditions for the asymptotic stability of synchronized inertial oscillator networks with delay. The corresponding delay master stability function (dMSF) is given in terms of the graph Laplacian spectrum, the delay τ , as well as dynamical parameters of the model. For delays caused by processing times, our results offer a complete analytic solution to the question of asymptotic stability. The dMSF comprises a finite number of easily evaluated critical conditions that hold for any network topology. Particularly, in an important case which covers our central application of renewable inverter-based power grids, the conditions further simplify to a single stability criterion involving just the maximum graph Laplacian eigenvalue.
Main application. We begin with an inverter-based power grid model to exemplify how we obtain a concise condition for a major application of oscillator networks. Due to the energy transition, power grids currently undergo substantial structural and dynamical changes, threatening stable synchronization of the AC voltage fre-quency [20,21]. Characterized by a large share of volatile distributed generation units, e.g. solar or wind power plants, future energy networks will require novel control approaches like grid-forming power inverters to maintain stability [22,23]. As this involves measurements and processing, delays are expected to play a critical role [24][25][26]. Understanding their influence on stability is thus vital to ensure security of supply and prevent blackouts.
Specifically, we consider frequency dynamics in a droop-controlled inverter grid [27] where the steady-state power flow between two nodes depends on the sine of their phase difference, Here ϕ i (t) denotes the phase angle of the i-th inverter (oscillator) and ∆ϕ τ ji .α > 0 andβ > 0 are the inertia-specific damping and droop constants, respectively; P d i represents the desired active power set points. Elements of the weighted adjacency matrix (K ij ) may be interpreted as the maximally transmittable power values along transmission lines in the network [28] (details in SI [29]).
We find that a synchronous state of Eq. (1) is asymptotically stable if and only if where y 1 ∈ (0, π]. This exact stability condition depends on network structure only via the largest eigenvalue λ N of the effective Laplacian matrix L (see below). Notably, it suffices to compute precisely one unique characteristic root y 1 of the linearized spectrum associated with Eq.
(1). We discuss this result further after deriving the general approach. Derivation. We consider a nonlinear dynamical system of N coupled oscillators on a network. All oscillators (nodes) have inertia, obeying a Newtonian law of motion. The state of the i-th oscillator at time t is given by the phase angle ϕ i (t) and angular frequency deviation ω i (t) ≡φ i (t) in a reference frame co-rotating with a coherent frequency Ω. Let the time evolution of the global state be governed by a set of second-order DDEs containing a discrete, constant delay τ > 0, for i ∈ {1, . . . , N }. Here time arguments are abbreviated as ϕ i ≡ ϕ i (t) and ϕ τ i ≡ ϕ i (t − τ ); furthermore ∆ϕ ij ≡ ϕ j − ϕ i and ∆ϕ τ ij ≡ ϕ τ j − ϕ τ i . The real scalar functions f i and f τ represent undelayed and delayed isolated dynamics, respectively. Unlike identical oscillators, f i may differ from node to node by an additional constant c i ∈ R, which accounts for heterogeneous driving forces. In the interaction term, g denotes undelayed coupling dynamics, whereas g τ describes interactions with a coupling processing delay. The strength of the coupling as well as the network topology are stored in the weighted adjacency matrix A ∈ R N ×N , with A ij > 0 if nodes i and j are connected and 0 otherwise. Here, we consider undirected graphs without self-loops.
The delay considered in Eq. (3) is a processing delay which arises, for example, in engineered systems with feedback control due to measurement and processing times. Contrarily, transmission or communication delays in diffusive coupling of the form g τ (x τ j − x i ) require separate treatment (see SI).
Traditional MSFs require complete synchronization, i.e. all oscillators move in phase with frequency Ω. Then, Jacobians evaluated on the synchronization manifold are identical for all nodes [13]. To achieve a block decomposition similar to MSF for phase synchronization, we assume: 1) The Jacobian matrices of all local functions f i and f τ , evaluated on Z, are identical.

2)
The Jacobians of all coupling functions g and g τ , evaluated on Z, are edge-independent except for a pre-factor w ij (∆ϕ * ij ) = w ji (∆ϕ * ji ) ∈ R which may depend on the fixed point ϕ * . We note that the following procedure also holds for more general coupling functions g(ϕ i ,φ i , ϕ j ,φ j ) and g τ (ϕ τ i ,φ τ i , ϕ τ j ,φ τ j ) if their first partial derivatives are antisymmetric with respect to the exchange of i and j. Nonetheless, we present the widely applied diffusive form here and refer to the SI for more information.
We define the effective Laplacian matrix L of the linearized network model such that This matrix is symmetric, positivesemidefinite, and consequently diagonalizable. In the spirit of MSF, we now transform coordinates into the space spanned by the eigenvectors of L, with corresponding eigenvalues λ k , k ∈ {1, . . . , N }. Diagonalization does not affect the Jacobians (which are node-independent by assumption after absorbing w ij in the adjacency matrix), such that the system of DDEs decomposes into N blocks given in terms of λ k , Here the set of θ k denotes (small) phase angle deviations from ϕ * expressed in the transformed coordinates. The coefficients are given by elements of the Jacobian matrices (see SI); in the following we suppress their dependence on λ k . The stability of a synchronous state ϕ * depends on the real parts of the roots of the characteristic equation associated with Eq. (4), The exponential polynomial H features an infinite number of complex roots; all must have negative real parts for asymptotic stability. We now assume that the delay τ appears either in the time argument of the phases ϕ i or of the frequencies ω i but not in both. For these cases, Bhatt and Hsu [19] derive necessary and sufficient stability conditions for scalar second-order DDEs, determined by a finite number of decisive roots within the infinite spectrum of H (see also [30,31]). After the block decomposition outlined above, we may transfer these conditions to inertial oscillator networks to obtain delay master stability conditions.
If only the coupling is delayed (f τ = 0), the longitudinal eigenvalue λ 1 = 0 describes dynamics within the synchronization manifold and asymptotic stability is determined by the N − 1 transversal directions [13]. Contrarily, if f τ = 0, all k must be considered for stability analysis. We define the transversal set N , which is {2, . . . , N } for f τ = 0 and {1, . . . , N } otherwise.
Substituting z = iy in Eq. (5), H separates into a real and an imaginary part. First, we consider the case of phase delay (a τ k = 0). Let a k > 0 and −a k < b k τ ≤ 0. Then, for each k, there exists one decisive root y 1,k ∈ (0, π] which solves the imaginary part of Eq. (5) [19]. A synchronous fixed point of Eq. (3) with phase delay is asymptotically stable if and only if, for all k ∈ N , with R k (y) : In the frequency delay case (b τ k = 0), we assume a k > 0 and b k > 0. Here we examine positive solutions y k of the real part of Eq. (5). The first positive root y k,0 lies in the interval (0, π/2), and one root y k,m is situated in each πinterval (mπ − π/2, mπ + π/2) for m = 1, 2, . . . . Of these roots, the decisive roots y * k and y * * k are found according to it is necessary and sufficient for asymptotic stability of a synchronous state that, for all k ∈ N , where R k (y) is defined beneath Eq. (6). Though finding decisive roots appears more complicated than in the previous case, it turns out that we must find at maximum 2N roots in total, within known intervals. Details are provided in the SI. Phase delay. Our motivating example introduced in Eq. (1) illustrates the case of phase delay. Here the root y 1,k is k-independent and we must only consider the largest eigenvalue λ N of L. We may write the resulting stability condition (Eq. (2)) as a dMSF σ, A combination of λ N and τ is stable if and only if σ < 0. Monotonicity arguments for y 1 (τ ) prove the existence of precisely one critical delay τ c = τ c (α,β, λ N ). A synchronous state that is stable without delay will remain asymptotically stable for all τ < τ c and is unstable for all τ ≥ τ c .
To visualize σ, we first choose a star topology as in Ref. [32]. A producer in the center with steady-state power production P + = 3P 0 is connected to three consumers with P − = −P 0 via transmission lines of equal capacity K 0 (Fig. 1a). Due to the symmetry of the configuration, we obtain three distinct Laplacian eigenvalues. The phase delay case depends on λ N only; thus we get a single curve with σ = 0 at τ c (Fig. 1b). For typical parameter values, τ c ≈ 40 ms is about twice the 50 Hz oscillation period.
Frequency delay. Stability in the presence of a frequency delay is qualitatively different. We show this by applying our approach to the decentral smart grid control (DSGC) scheme [32,33]. The model incorporates electricity price dynamics by locally relating the price to the current grid frequency, motivating producers/consumers to adapt their feed-in/consumption to the currently available power supply. The continuous measurements required for this smart grid regulation induce a local processing delay at each node. The model reads  [32]φ where P i is the produced/consumed power at node i and K = (K ij ) denotes the weighted adjacency matrix as before. Next to the damping α > 0, the price elasticity γ > 0 acts as a second, delayed damping term. For the DSGC model we obtain N delay master stability conditions Here it is not a priori identifiable which root y * k = y * k (λ k , τ ) determines stability for a given τ ; we must regard all Laplacian eigenvalues λ k . The root y * * k is not relevant because γτ > 0. Analogous to the previous case, we may formulate Eq. (10) as a dMSF σ(λ k , τ ). Then, the system is stable in all regions where σ max (τ ) := max Calculating σ max for the star topology (Fig. 1a), we now have contributions from all three distinct eigenvalues λ k as depicted in Fig. 1c. In addition to a stable regime beginning at τ = 0, there exist further windows of stability for larger delays, corroborating prior results based on numerical analysis [32].
Networks. For frequency delays, we conjecture that the length of the first stability window, extending from τ = 0 to a critical delay τ c , is determined by λ N . In the limit of an infinite, heterogenous graph, its Laplacian spectrum may be expected to become quasi-continuous, such that further stability windows vanish. Thus, in both frequency and phase delay, the maximum eigenvalue of the effective graph Laplacian plays a crucial role. Several bounds and estimates in terms of network characteristics have been published for the largest Laplacian eigenvalue of a weighted graph, e.g. [34]. Therefore, our method may provide insight even without explicitly diagonalizing L, which could be practical particularly for large systems.
Finally, we explore how the critical delay τ c (end of first stability window) depends on the network structure in both models. As a versatile example, we generate Watts-Strogatz networks [35] with different rewiring probabilities p, mean degrees k, and number of nodes N (e.g. Fig.  2a). Consumers and producers are placed alternately on the original ring graph. The results for varying p, k, and N (Figs. 2b-d) show that the critical delay decreases with an increasing number of nodes and edges as well as randomness of the system, which emphasizes the significance of delays for the design and control of real-world dynamical systems. However, the decline is smaller for the DSGC model compared to the droop-controlled inverter model.
Conclusion. In this Letter, we present an analytical approach to assess the asymptotic stability of synchronized inertial oscillator networks with delayed dynamics. Specifically, we consider processing delays either in the phase or in the frequency. We show how to extend the master stability formalism to an arbitrary lag time, obtaining dMSFs in terms of the delay τ and eigenvalues of the effective graph Laplacian L. Unlike MSF, we more generally consider phase synchronization and allow for constant inhomogeneities in the local dynamics of otherwise identical oscillators. A block decomposition of the linearized model yields necessary and sufficient stability conditions which deliver an analytic expression for the dependence on network structure. These criteria involve maximally 2N decisive roots. In contrast, previous numerical asymptotic stability analyses rely on randomly computing a significantly larger number of characteristic roots without being certain that all decisive roots have been found. Illustrating our approach, we consider two concrete models for renewable power grids as our main applications. Notably, we are able to boil down the problem of stability to a single condition in the case of the droop-controlled inverter model. We therefore believe that our method could contribute to the development of design criteria for future energy systems. Generally, our results advance the stability analysis of dynamical systems on complex oscillator networks complying with Eq. (3).
Due to its increasing importance in real-world applications, delay stability in complex systems remains an important topic with many open challenges. Our work opens a new analytic approach on this subject. From the perspective of power grids, it is crucial to also tackle non-identical oscillators, to consider multiple delays, and combined phase and frequency delays. In the wider context of physical systems, it is also highly interesting to study more general stability criteria for diffusive coupling with transmission delays.
The authors acknowledge the support of BMBF, Con-dynet2 FK. 03EK3055A. This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -KU 837/39-1 / RA 516/13-1. All authors gratefully acknowledge the European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research and the Land Brandenburg for supporting this project by providing resources on the high performance computer system at the Potsdam Institute for Climate Impact Research.

I. COEFFICIENTS OF SECOND-ORDER DDE BLOCKS
The coefficients in Eq. (4) of the main text are given by Elements of the delayed Jacobians are written analogously as F τ ϕ etc.

II. GENERAL DERIVATION
In the Letter, we outline the derivation of our approach based on the inertial oscillator model described by Eq. (3) of the main text. Here we present a more elaborate variant (discussed in detail in Ref. [36]). This allows us to illuminate the underlying assumptions of our method and discuss why including communication delays involves a strong restriction.
Recall that we formulate the oscillators' dynamics in a reference frame co-rotating with the frequency Ω. A synchronous state where all oscillators are entrained to frequency Ω corresponds to a fixed point characterized by Instead of the second-order form stated in Eq. (3) of the main text, we may equivalently express our inertial oscillator network model as a set of first-order DDEs by treating the phase angles ϕ i (t) and angular frequency deviations ω i (t) =φ i (t) as two independent variables for each node. We thus define the vector x i ≡ (ϕ i , ω i ) and writeẋ Here we denote vector-valued functions (R 2 → R 2 ) by bold letters, while f i , f τ , g, and g τ will remain the scalar functions introduced in the main text; e.g. f i (x i ) = (ω i , f i (ϕ i , ω i )) . In contrast to Eq. (3), this form now includes a communication delay via the function g 0τ (x i (t), x j (t− τ )). Furthermore, note that the coupling functions g 00 , g τ τ , and g 0τ may depend on x i and x j in an arbitrary fashion.
A communication delay of the type above may describe transmission or propagation lags between nodes. The intuition is that a change of node i at time t depends on the history of connected nodes because it takes the time τ until a signal from a node j reaches node i.
We now linearize Eq. (S.3) around a fixed point x * = (x * 1 , . . . , x * N ) defined by the conditionsẋ * i = (0, 0) for all i. The set of fixed points constitutes the phase synchronization manifold Z, where Jacobian matrices, all evaluated at the fixed point, are written in short notation, . Due to the relation ω i =φ i between coordinates of the vector x i , some elements of the Jacobians are immediately zero or one. Particularly, where ∂ m f denotes the partial derivative of the function f by the argument m, evaluated at the fixed point. In the same manner, We emphasize that these Jacobians depend on the fixed point. For the phase synchronization manifold, this implies that the Jacobians may differ for different i and j.

A. Antisymmetric coupling
Assume now that 1. there is no communication delay (i.e. g 0τ = 0) and 2. the linearized coupling between two nodes i and j is antisymmetric, that is, This is fulfilled by the model discussed in the main text (Eq. (3)), where we have diffusive coupling. The antisymmetry requirement will allow us to write the problem in terms of the effective graph Laplacian matrix L.
Our goal is to decouple local information about the dynamics of single nodes from global terms characterizing the network as a whole. Mathematically, this is achieved when local 2 × 2 matrices and global N × N matrices factorize into a Kronecker product (symbolized by ⊗).
If we have complete synchronization (i.e. all nodes oscillate with identical frequency and phase angle), the Jacobians are homogeneous for all i, j, resulting in immediate Kronecker factorization. This is not true for the more general case of phase synchronization. To achieve a decomposition nonetheless, we impose the following restrictions in analogy to the main text: 3. The Jacobians of the local functions f and f τ , respectively, evaluated on the phase synchronization manifold, are identical for all nodes: 4. The Jacobians of the coupling functions g 00 and g τ τ , respectively, evaluated on the phase synchronization manifold, are identical for all i, j up to a prefactor w ij (x * i , x * j ) ∈ R which contains all dependencies on the fixed point. It has the property w ij = w ji > 0; With these assumptions, all dependencies on the fixed point may be absorbed in the effective adjacency matrix A with entries A ij := w ij A ij . Since we assume antisymmetric coupling (assumption 2), we may furthermore replace A by the the effective graph Laplacian matrix L given by L ij := δ ij j A ij − A ij . Now, the set of linearized DDEs readṡ or, in vector notation for the entire system, η = (η 1 , . . . , η N ) , Here I N is the N -dimensional unit matrix. According to assumption 4, L is symmetric and therefore diagonalizable. Switching to a basis B of eigenvectors via the coordinate transform ξ = [T B ⊗ I 2 ]η, we diagonalize L = T −1 B ΛT B to obtain the diagonal matrix of its eigenvalues, Λ = diag(λ 1 , . . . , λ N ). This leads to a block-diagonal form; each two-dimensional block is given by the equatioṅ where k ∈ 1, . . . , N . (Note that we have switched index from i to k to emphasize that ξ represents the state vector in the eigenbasis B.) Similar to the master stability formalism, we have thus decomposed the problem into blocks which vary only in the effective Laplacian eigenvalue λ k . All eigenvalues are non-negative and λ 1 = 0 because, according to the properties of an undirected graph's Laplacian matrix, L is positive-semidefinite.
Recalling that the two components of the vector x i are related via ω i =φ i , we may write Eq. (S.9) in second-order form. The vector φ of all phase angles is given by the linear combination φ = k v k θ k , where v k is an eigenvector of L. In terms of the phase angles in the transformed coordinates, θ k , we obtain a second-order DDE, The coefficients are given in Eq. (S.1) as functions of λ k and the Jacobians F , F τ , G, and G τ τ . Following Bhatt and Hsu [19], we state stability criteria in terms of these coefficients, as discussed in the main text.

B. Communication delays
The derivation above relies on the assumption of antisymmetric coupling between connected nodes i and j (assumptions 1 and 2). This permits rewriting the problem in terms of the effective Laplacian matrix L, which leads to a decomposition in its eigenvalues.
Let us briefly consider symmetric or asymmetric coupling which does not satisfy assumptions 1 and 2. Notably, a communication delay immediately destroys the antisymmetry because the coupling function g 0τ (x i , x τ j ) evaluates its arguments at different times. In this case, the problem cannot be expressed in graph Laplacian form. Instead, we consider an effective adjacency matrix A and an effective degree matrix D. These matrices must then commute to allow a block decomposition.
In the following, let g represent any of the functions g 00 , g τ τ , or g 0τ . When linearizing the inertial oscillator model, we obtain one Jacobian D 2 ij g(x * i , x * j ) containing derivatives with respect to the j-th node, and a second Jacobian D 1 ij g(x * i , x * j ) with derivatives by x i . Replacing assumption 4, we assume: 5. The adjacency matrix A and the Jacobian D 2 ij g(x * i , x * j ), evaluated on the phase synchronization manifold Z, factorize into the direct product A ⊗ G (2) of an effective N × N adjacency matrix A and a universal 2 × 2 Jacobian G (2) , such that the local matrix G (2) is the same for all nodes and only A depends on the indices i, j of the network. Similarly, A and D 1 ij g(x * i , x * j ) factorize into the direct product A ⊗ G (1) . We require that A and A are symmetric matrices.
Since A and A are symmetric by assumption, they are diagonalizable. We define the effective degree matrix D, Now we assume: 6. The matrices A and D commute, i.e.
This implies that they are simultaneously diagonalizable. In that case, the DDE (Eq. (S.4)) decomposes into blocks in terms of the eigenvalues of A and D. Then, the coefficients in Eq. (S.10) are functions of these eigenvalues instead of the eigenvalues λ k of the graph Laplacian matrix L, and our method can be applied in analogy to the antisymmetric case. Regular graphs with homogeneous weights present a special case where [A, D] = 0. For more complex network topologies, however, assumption 6 is generally not satisfied.

III. DECISIVE ROOTS
In the case of phase delay, there is in principle one uniquely defined decisive root for each k. Thus, the stability analysis is generally based on calculating maximally N decisive roots for a system of size N . In contrast, the frequency delay case generally requires identification of two decisive roots y * k , y * * k for each k. The value of these roots depends on the delay and on the coefficient b k .
Assume the value √ b k τ 2 is located in the interval (jπ − π/2, jπ + π/2), where j is a positive integer. If j is an odd number, the closest decisive root y * k is in the same interval. We then find the root y * * k either in the π-interval to the right (j + 1) or to the left (j − 1). To find all decisive roots for a given k, we must therefore calculate three roots and, among them, compare the two possible candidates for y * * k . If j is an even number, then y * * k is located in the same interval and y * k is one of the two roots found in each of the adjacent intervals. In any case, all decisive roots must lie within a distance of 3π/2 from the value √ b k τ 2 (see figure S.1).
In conclusion, for frequency delays, we must calculate at maximum 3N characteristic roots to find in total 2N decisive roots that the stability conditions are based on. In the DSGC model presented in the main text, the root y * * k turns out to be irrelevant; thus it suffices to calculate at maximum 2N characteristic roots.

IV. DROOP-CONTROLLED INVERTER MODEL
In this section, we provide further detail on the renewable inverter-based power grid model with processing delay, considered in the main text as a central application of the phase delay case. For additional information on the theoretical study of power systems we refer to [28,[37][38][39].
To model a renewable power system with phase delay, we describe the dynamics of grid-forming inverters (represented by nodes of the network) using the swing equation [28,37], Here ϕ(t) denotes the phase angle, m i the inertia, a i the damping constant, and P d i the desired power set point of the i-th node. The set point is positive for production and negative for consumption. Furthermore, we write the electrical power at node i as where U 0 denotes the AC voltage amplitude which is assumed constant throughout the system, i.e. U 0 = U i ∀i [39]. B ij represents the susceptance of the transmission line between nodes i and j (we may choose its value to be zero if i and j are not directly connected). Eq. (S.12) expresses a common choice in the literature to model steady-state power flow [28]. It follows when neglecting losses (purely inductive power lines) and assuming that all phase differences |ϕ j − ϕ i | < π/2. According to Eq. (S.12), the power flow along a transmission line between two nodes depends on the phase angle difference between them. Let us suppose that the state of this power line (an edge in the network) enters the frequency control, yet with a processing delay. One way to model this would be delayed coupling, where the phase difference is evaluated at time (t − τ ).