Fluctuation response patterns of network dynamics – An introduction

Networked dynamical systems, i.e., systems of dynamical units coupled via nontrivial interaction topologies, constitute models of broad classes of complex systems, ranging from gene regulatory and metabolic circuits in our cells to pandemics spreading across continents. Most of such systems are driven by irregular and distributed fluctuating input signals from the environment. Yet how networked dynamical systems collectively respond to such fluctuations depends on the location and type of driving signal, the interaction topology and several other factors and remains largely unknown to date. As a key example, modern electric power grids are undergoing a rapid and systematic transformation towards more sustainable systems, signified by high penetrations of renewable energy sources. These in turn introduce significant fluctuations in power input and thereby pose immediate challenges to the stable operation of power grid systems. How power grid systems dynamically respond to fluctuating power feed-in as well as other temporal changes is critical for ensuring a reliable operation of power grids yet not well understood. In this work, we systematically introduce a linear response theory (LRT) for fluctuation-driven networked dynamical systems. The derivations presented not only provide approximate analytical descriptions of the dynamical responses of networks, but more importantly, also allow to extract key qualitative features about spatio-temporally distributed response patterns. Specifically, we provide a general formulation of a LRT for perturbed networked dynamical systems, explicate how dynamic network response patterns arise from the solution of the linearised response dynamics, and emphasise the role of LRT in predicting and comprehending power grid responses on different temporal and spatial scales and to various types of disturbances. Understanding such patterns from a general, mathematical perspective enables to estimate network responses quickly and intuitively, and to develop guiding principles for, e.g., power grid operation, control and design.


Introduction
Networked dynamical systems abound around and in us.From the circuits supporting metabolism and gene regulation in our cells to the neural networks in our brains, from the power grids that supply electric energy to most of our technical infrastructure to the internet that connects our computers, all of these systems are driven externally, often by irregular, time-dependent and spatially heterogeneous signals.How networked dynamical systems respond to such perturbations, driving signals, or other types of time-dependent inputs is hardly understood to date.In this article, we offer a general introduction to the basic theory of analyzing response features of networked dynamical systems by linear response theory and focus on applications in the realm of power grids.
A reliable supply with electric power fundamentally underlies most aspects of modern society.As the shares of renewable electric energy supply grow and consumer dynamics is increasingly influenced by digital technologies, fluctuations impinge on power grids, making them intrinsically driven, non-equilibrium systems, with distributed and often non-stationary response dynamics Witthaut et al. (2022).If fluctuations become large, they may destabilize grid dynamics, affect normal operations of parts of the grid, cause cascades of failures or even total blackouts [UCTE (2007), Wilde (2020), ENTSO-E (2021)].To predict, control or mitigate the fluctuating and distributed responses of such networks resulting from input (and output) power fluctuations, we need to understand their nature in network dynamical systems.
How can we systematically characterize fluctuating responses that are distributed across meshed networks?How can we predict their influences and at which nodes in a network is their impact most profound?Linear response theory (LRT) relates sufficiently small time-dependent driving signals to time-dependent responses.The theory approximates a resulting system dynamics near some operating point to first order in the strength of the driving signals.It has traditionally found applications across physics, chemistry and engineering [Kubo (1957), Cammi and Mennucci (1999), Ikeguchi et al. (2005), Ruelle (2009), Majda and Wang (2010), Pan et al. (2020)], and recently also in power grid models [Dörfler et al. (2013), Witthaut et al. (2016), Kettemann (2016), Manik et al. (2017), Tamrakar et al. (2018), Haehne et al. (2019), Tyloo et al. (2019)].However, a framework of linear response theory uncovering spatiotemporal response patterns in systems with multi-dimensional dynamics of units that simultaneously interact via intricate network topologies is not yet well established.
In this article, we introduce a general formalism of linear response theory for networked dynamical systems and demonstrate its applications in stationary and non-stationary models of power grids.Specifically, in Sec. 2, we present the main ideas of the linear response theory for networks by first analyzing the linear responses of networked dynamical systems with the most general settings (Sec.2.1) and then boiling down step-by-step to a specific LRT which provides a direct link between the temporal and the spatial features of the dynamic network responses (Sec.2.2).In Sec. 3, we demonstrate the LRT for networks by applying it to two models of power grids: a stationary model for the DC approximation of the AC power flow distributions in power grids (Sec.3.1), and a non-stationary model, the (second-order) oscillator model, commonly used to describe the dynamics of the high-voltage AC power transmission networks (Sec.3.2).The next two sections, Sec. 4 and Sec. 5, focus on the extraction and the interpretation of spatiotemporal response patterns of power grids from the LRT, as well as the analytical techniques needed therein.In Sec. 4, we elaborate how distinctive steady-state response patterns emerge in three frequency regimes (Sec.4.1), how the transient spreading pattern of a perturbation entangles with the underlying network topology (Sec.4.2), and how LRT helps to estimate the long-term risk of individual nodes from external fluctuations (Sec.4.3).In Sec. 5 we summarize the role of LRT in uncovering such patterns by comparing the analyses for patterns in the steady-state response vs. the transient responses (Sec.5.1), the patterns in the deterministic responses to a given perturbation vs. the cumulative responses to a random signal (Sec.5.2), and the small responses close to the base operation state vs. the large responses further away (Sec.5.3).In the last section (Sec.6), we point out several directions for extending the present theory for a better understanding of the response dynamics of networked dynamical systems and for a safer and more reliable operation of future power grids.

Main ideas of the linear response theory for networks
In this section we introduce the main ideas of the linear response theory framework.Sec.2.1 formulates the LRT on a general model of networks of one-dimensional dynamical units.Sec.2.2 highlights the linear operators, represented as matrices, arising in the LRT for specific network settings and how basic information of the underlying network system such as its topology and features defining the dynamically determined operating point enter those operators.In Sec.2.2 we also derive the explicit linear responses of the general network model of one-dimensional units introduced in Sec.2.1 and extend the results for networks of higher-order nodal dynamics.

General formulation of LRT for networked dynamical systems
We now illustrate the main idea underlying LRT by starting with a general setting of networked dynamical systems where each unit's dynamics is one-dimensional.We consider a dynamical process involving N interacting variables, whose state is represented by a vector x = (x 1 , • • • , x N ) ∈ R N .The autonomous dynamics of the N -dimensional dynamical system is governed by where f : R N → R N is a function that in general depends on the states x of all units and does not explicitly depend on time.Let us consider a system that exhibits a fixed point at x = x * where f (x * ) = 0, and is influenced by an external dynamic driving vector D(t) ∈ R N at the fixed point x * .We investigate how the system dynamically responds to D(t), i.e. how the system's deviation X(t) := x(t) − x * from the fixed point evolves in time.In general, the dynamics of the system's response X follows where both functions f and D can be highly nonlinear so that an exact analytical solution of the system's response X typically does not exist or is unknown.Important information about the response dynamics (2.2) is given by the stability operator at the fixed point obtained from the linearization of the function f at x = x * , i.e. the Jacobian matrix J with J ij := ∂fi ∂xj | x=x * .The system's response dynamics thus approximately follows the linearized differential equation Ẋ = J X + D(t). (2. 3) The signs of the real parts of the Jacobian eigenvalues w [ℓ] (with ℓ being the index) indicate whether the fixed point x * is linearly stable (Re[w [ℓ] ] < 0), unstable (Re[w [ℓ] ] > 0) or neutrally stable (Re[w [ℓ] ] = 0) in the eigendirection/eigenspace corresponding to the eigenvalue w [ℓ] .Below, we typically focus on the dynamics near stable fixed points, where all Re[w [ℓ] ] < 0 (or marginally stable ones where one Re[w [ℓ] ] = 0 if J is a Laplacian operator) to justify the assumption that the solution X(t) of the linearized equation (2.3) reasonably well approximates the full solution of (2.2).
Let us now focus on systems with pairwise interactions between the N units such that the function f i (x) controling the time evolution of unit i can be written as (2.4) where h i denotes the intrinsic dynamics depending on the variable x i itself and the coupling term K ij g ij (x i , x j ) represents the pairwise interaction between variable x i and x j with i ̸ = j.Here K ij ∈ R + 0 is the non-negative coupling strength and g ij is the coupling function depending on the state of the pair of variables (x i , x j ) and g ij (x i , x j ) ̸ ≡ 0.
Remark 2.1 Pairwise interactions (2.4) induce an interaction structure of the dynamical system (2.1) that can be represented by a graph G(V, E) where the set of vertices V consists of N variables x 1 , • • • , x N and the set of edges E consists of all node pairs with the pairwise coupling strength being nonzero, i.e.E = {(i, j) The linearization of the pairwise interaction (2.4) yields the Jacobian matrix of the response dynamics (2.2) of the driven networked system (2.5 b) It is clear that the topology of the interaction network G explicitly enters the Jacobian matrix through the matrix of the coupling strength K ∈ R N ×N with its ij-th element being defined as K ij , which is effectively a weighted adjacency matrix of graph G.
The combination of the linearized dynamics of a general N -dimensional networked dynamical system (2.1) near a fixed point (2.3) and the arising Jacobian matrix that explicitly depends on the system's topology (2.5) due to pairwise interactions (2.4) provides a general form of LRT for networked dynamical systems.For a specific system with given forms of or constraints for the intrinsic nodal dynamics h i (x i ), coupling strengths K ij and coupling function g ij (x i , x j ), the solution of the matrix equation (2.3), X(t), can be characterized by evaluating the specific spectral properties of the Jacobian matrix J .Thereby the LRT provides a powerful tool to describe the dynamic response of a networked system temporally and spatially at once: the solution X(t) explicitly depends on time, and at the same time its relation to the topology-dependent Jacobian matrix can be exploited to reveal the temporally and spatially distributed patterns of the network response.

Arising linear operators and network topology
As discussed in the last section, the dependence of the Jacobian matrix of network response dynamics on a weighted adjacency matrix of the underlying interaction network gives the first hint on how temporal and spatial features of dynamic network responses are intertwined.In this section we further show that, under a few commonly satisfied conditions such as diffusive coupling between units, interesting results of network response dynamics emerge.Especially, another important graph-theoretical matrix, the Laplacian matrix, arises in the network response dynamics (2.3), providing us powerful tools for characterizing the spatiotemporal patterns in dynamic network responses.
Diffusive coupling is a very common type of coupling present in many physiological and chemical systems [Hale (1997), Larter et al. (1999), Postnov et al. (2006), Casagrande (2006), Stankovski et al. (2017)], in particular also appearing in the Kuramoto model [Kuramoto (1984)] and its variations [Filatrella et al. (2008), Acebrón et al. (2005)].A diffusive coupling function gij mediating the interaction between a pair of nodes (i, j) is characterized by its dependence on the state difference x j − x i of the node pair, i.e. g ij (x i , x j ) = gij (x j − x i ) in (2.5).For notational simplicity, we again denote the functions gij just by g ij .
Proposition 2.1 (Stability of diffusively coupled networks: a special case) A networked dynamical system with evolution function (2.4) and interaction network G(V, E) is at least neutrally stable at its fixed point x = x * if a) the intrinsic nodal dynamics Proof The diffusive form of the coupling function g ij (x j − x i ) yields a useful relation . With this particular symmetry, the Jacobian matrix (2.5) of the system at the fixed point x = x * takes the following form where (2.7)According to the Gershgorin circle theorem, every eigenvalue of the Jacobian matrix lies within at least one of the N Gershgorin discs in the complex plane.Relation (2.7) indicates that all N Gershgorin discs lie in the left half of the complex plane, i.e. in {z ∈ C | Re (z) ⩽ 0 }, because the centers of the discs (J ii , 0) lie on the negative real axis and the radius of the discs r i = k̸ =i |J ki | ⩽ |J ii |.The discs touch the imaginary axis from the half plane {z ∈ C | Re (z) ⩽ 0 } only if β i = 0. Therefore, all Jacobian eigenvalues can only have non-positive real parts and consequently the networked dynamical system is at least neutrally stable at the fixed point x = x * .Remark 2.2 (Homogeneous nodal dynamics and graph Laplacian) We remark that in case of identical intrinsic nodal dynamics at all nodes, a weighted Laplacian matrix L of the interaction graph explicitly enters the linearized response dynamics of the network (2.3).Assuming β i = β ∈ R for all nodes i, we can express the Jacobian matrix as where the weighted graph Laplacian L is defined as (2.9) Here K ij γ ij is considered a weight of edge (i, j), containing the coupling strength K ij and the linearized coupling function γ ij at the fixed point.
Remark 2.3 (Symmetry and linear responses in Laplacian eigenbasis) In general, the interaction network can be directed, meaning that for an edge (i, j) the coupling strength K ij and the derivative of the coupling function γ ij at the fixed point, i.e. the sensitivity of the diffusive coupling function g ij (x j − x i ) to a change in the state difference x j − x i , can differ from their counterparts K ji and γ ji for the edge (j, i) with the opposite direction.This asymmetry leads to an asymmetric weighted graph Laplacian (2.9).However, for undirected networks with symmetric strengths (K ij = K ji ) and symmetric sensitivities of coupling functions (γ ij = γ ji ), or more generally, a symmetric combination K ij γ ij = K ji γ ji , the weighted graph Laplacian L is symmetric.Its eigenvectors thus form an orthogonal basis which allows us to solve for the linear network responses in (2.3) and (2.8) by expressing the response vector X(t) in terms of the eigenvalues and eigenvectors of the Laplacian.
We first analyze a system that is perturbed by only one sinusoidal signal with magnitude ε > 0 and frequency ω > 0 at node k, i.e.D (k) ωt+φ) , where δ ik is the Kronecker delta with δ ik = 1 if i = k and δ ik = 0 otherwise.We solve for the linear network response vector (2.10) Expressing the response vector in the constant eigenbasis of Laplacian (2.11) and the time-dependent coefficients c [ℓ] (t) and exploiting orthogonality of the Laplacian eigenvectors, we obtain the ordinary differential equations ωt+φ)  (2.12) for the coefficients.Here the N Laplacian eigenvalues λ [ℓ] and eigenvectors v [ℓ] are indexed . If the graph is connected, the zero eigenvalue is unique, such that all other eigenvalues are positive real numbers [Bapat (2010)].The solution of differential equation (2.12) for the coefficients is (2.13)where X (k) 0 denotes the initial response vector at t = 0 given node k is perturbed.The linear network response is thus given by (2.11) and (2.13).
Remark 2.4 (LRT for distributed arbitrary perturbations) In case perturbations with arbitrary temporal structures are distributed across the network, that is, multiple elements of the perturbation vector D(t) are arbitrary time series, the linear network response can be obtained by summing up the single-signal single-frequency response (2.11 and 2.13) over all frequency components ω and all perturbation signals k by resorting to the linearity of the dynamics (2.3) (2.14) Here X (k) ω=0 (t) denotes the linear response to a constant (ω = 0) perturbation ε, which shares the same form as X (k) (ω, t) with ω = 0 for β > 0 but not for β = 0.Moreover, the sum becomes an integral for continuously distributed frequencies, with X (k) (ω, t) replaced by the associated response density per unit frequency.
For systems where nodal damping vanishes (β = 0), the coefficient for the 0-th eigenmode c [0] diverges for a constant perturbation as ω → 0, so it needs to be solved separately and takes the form of k t, inducing a linear drift and thus unbounded growth with time.As the dynamics LRT describes is only approximate for nonlinear dynamical system, unbounded growth typically induces the approximation of the linear response to the real system to break down, i.e. become useless in practice due to larger errors between approximate and exact solution.
Remark 2.5 (LRT for higher-order nodal dynamics) In the above paragraphs we discussed the main ideas of the LRT for networked dynamical systems with first-order nodal dynamics.For more general systems with second-or higher-order nodal dynamics, the straightforward relation (2.8) between the Jacobian matrix and a weighted graph Laplacian does not hold any more.Nevertheless, a symmetric weighted graph Laplacian still arises in the linearized response dynamics for diffusively-coupled undirected networks with symmetric coupling strengths and symmetric sensitivities of coupling functions as discussed in Remark 2.3.If the higher-order time derivatives of the state variables has homogeneous coefficients for individual nodes, i.e. the response dynamics to a perturbation D(t) has the form of an explicit solution of the linear responses in the eigenbasis of L can still be obtained following the routine in Remark 2.3, if the corresponding ODEs for the time-dependent coefficients (2.16) are solvable.Here κ d ∈ R are constant coefficients and we use Euler's notation for derivatives, where D d t x denotes the d-th time derivative of variable x.
In summary, if a networked dynamical system consisting of N units (1) has a fixed point at x = x * , (2) in the neighbourhood of x = x * the intrinsic nodal dynamics h i (x i ) gives homogeneous non-positive feedback to the respective nodes, i.e.

and
(3) the coupling function g ij (x j − x i ) is diffusive and the coupling term's sensitivity to small changes at x = x * is symmetric and non-negative, i.e.
= γ ji ⩾ 0 for all edges (i, j), then i) the dynamical system is at least neutrally stable (Proposition 2.1) and ii) a symmetric weighted graph Laplacian (2.9) arises in the network response dynamics (2.10) and enables the expression of linear network responses in the Laplacian eigenbasis (2.11) and (2.13) (Remark 2.2 and Remark 2.3).As we will show in Section 4, the explicit dependence of the linear network response on Laplacian eigenvalues and eigenvectors provides a powerful tool to reveal how the dynamic responses spatially distributed across the network.

LRT for power grid models
We now discuss how the LRT for general networked dynamical systems introduced in Section 2 applies to stationary and non-stationary models of power grids and helps reveal static and dynamic responses of power grid systems to external perturbations.

LRT for the DC power flow model
We first demonstrate how LRT works in a minimal model, the DC power flow model, and how it helps to compute the systemic stationary response of a power transmission network to perturbations in power injections and withdrawals.
For common AC power grids, the full power flow analysis poses several challenges such as possible difficulties in finding a solution in ill-conditioned cases and the existence of multiple solutions due to the inherent nonlinearities [Milano (2010)].By linearizing the AC power flow equations at an operation point, the DC power flow model1 provides a relatively simple and computational inexpensive way to compute the power flows.
In an AC power transmission grid, the total complex power flow S j from unit j to unit i reads where U j = |U j |e ıθj and I ij = (U j −U i )/Z ij denote the voltage at node j and the current between nodes j and i, respectively.Both are expressed as complex numbers to reflect the oscillating nature of AC power generation.The asterisks in, e.g., I * ij indicate the complex conjugates (e.g., of I ij ).Moreover, Z ij = R ij + ıX ij denotes the impedance of the transmission line (i, j) between unit i and j, with R ij denoting the resistance and X ij the reactance of (i, j).The complex power flow S ij = P ij + ıQ ij consists of two parts, the active power P ij that results in net energy transfer and the reactive power Q ij that returns to the source in each cycle, not doing any work, but supporting the voltage stability of the power system [Machowski et al. (2008)].
The DC power flow model is based on the following assumptions on the parameters and the operating state of the power grid systems: Assumption 3.1 (Perfect voltage support) The voltage amplitude is constant and identical for each node in the power grid network, |U i | ≡ U for all i, and the management of the reactive power is not considered.Assumption 3.2 (Lossless lines) Transmission losses on the lines are neglected, implying that line resistances are negligible compared to line reactances: R ij /X ij → 0 for all lines (i, j).Assumption 3.3 (Low line loads) Loads on all transmission lines are low, that is, the voltage angle differences between all neighboring nodes are much smaller in magnitude than π/2 such that sin(θ j − θ i ) ≈ θ j − θ i and cos(θ j − θ i ) ≈ 1.
With the above mentioned assumptions in mind, the complex power flow (3.1) between neighboring nodes simplifies to Here the complex power flow S ij naturally reduces to the active power flow P ij since the imaginary part vanishes.The equation (3.2) resembles the expression of the direct current carried by a "resistor" X ij /U 2 influenced by a "voltage drop" θ j − θ i according to Ohm's law, hence the name "DC power flow model".
For an AC power grid system consisting of N units, the active power flow P i injected at unit i is the sum over all connected units.Equation (3.3) is the core of the DC power flow model as it yields the pattern of power flows K ij θ j across the grid network.Here we follow the notation introduced in Sec.2.1 and define the coupling strength as K ij = U 2 /X ij if there exists a transmission line between unit i and j and K ij = 0 if there is not.Denoting the nodal active power injections and nodal voltage angles as vectors, P := (P 1 , • • • , P N ) and θ := (θ 1 , • • • , θ N ), we express the linear relation between them by a weighted graph Laplacian L introduced in Section 2.2, i.e.
Here L is defined similarly as in (2.9),only with γ ij ≡ 1 for all edges (i, j).
Assuming that the power transmission network runs at a normal operation state where the voltage angles θ * are stationary at all nodes, the fixed voltage angle differences determine a specific power flow pattern P * across the network through the linear operator L such that P * = Lθ * .If the nodal power injections are perturbed as P (t) = P * + D(t) by a shift vector D(t) that in general is time dependent and reflects an increase or decrease of power generation for consumption at some of the nodes, the nodal voltage angles θ = θ * + Θ change accordingly through L due to the linear operation in (3.4).The response vector of the voltage angles Θ is given by where L + denotes the Moore-Penrose inverse of the weighted graph Laplacian L.
Remark 3.1 The weighted graph Laplacian L defined in (2.9) is singular since it always has an eigenvalue λ 0 = 0 with the corresponding eigenvector v 0 = (1, • • • , 1) satisfying Lv 0 = 0 by construction.Therefore, to compute the voltage angle shifts (3.5) we need the generalized inverse matrix L + , which can be computed by, e.g., the singular value decomposition.Alternatively, we can remove one dimension of the system by treating the phase θ k of one of the units k as the reference for voltage angle, i.e. by setting θ i → θ i − θ k for all i.If the network is connected, then the (N − 1)-dimensional Laplacian matrix is invertible and both of the matrix equations, (3.4) and (3.5), have a unique solution respectively.

LRT for the oscillator model of AC power grid dynamics
In this section we discuss how LRT applies for the oscillator model of AC power grids, a widely used model for analysing the dynamics of AC power grids, and thereby provides a way to accurately determine the high-dimensional dynamic responses of an arbitrary power grid network to fluctuating power injections and withdrawals.
The dynamics of the high-voltage AC power transmission networks is essentially captured by an oscillator model (or second order model) of AC power grids, of which synchronization in terms of networked dynamical systems have been initially studied in references Filatrella et al. (2008), Rohden et al. (2012) and Motter et al. (2013).This model allows for analytical understanding of the dynamics of power grids and has yielded fruitful research results over the past decade [Rohden et al. (2012), Motter et al. (2013), Dörfler et al. (2013), Witthaut et al. (2016), Tyloo et al. (2019)].As the name suggests, in the oscillator model, each unit of AC power grids, a synchronous machine, is represented by an oscillator and the power transmission lines are represented by the pairwise couplings between the oscillators.The normal operation state of a power grid corresponds to the synchronization of all oscillators, where all units rotate at the same frequency Ω m 0 corresponding to the nominal grid frequency Ω 0 = 2π × 50 or 2π × 60 Hz.For each unit in the oscillator model, a synchronous machine, any change of the angular velocity of rotation results from the imbalance of the torques acting on the rotor operated at the nominal grid frequency.Its dynamics is governed by the so-called swing equation [Kundur et al. (1994), Machowski et al. (2008)]: (3.6)where θ m denotes the mechanical rotor angle deviation from the rotating reference frame Ω 0 t, I denotes the moment of inertia of the rotor and the connected turbine, D m denotes the coefficient of the damping torque resulting from the velocity-dependent friction at the air gap between the rotor and the stator in the synchronous machine.T m and T e denote respectively the net mechanical torque and the counteracting electromagnetic torque acting on the rotor.For deviations of the angular velocity θm , the local frequency deviation is small compared to the nominal grid frequency Ω m 0 , i.e. ( θm +Ω m 0 ) −1 ≈ (Ω m 0 ) −1 , so that a torque T acting on the rotor can be expressed in terms of the power P of the synchronous machine as T = ( θm + Ω m 0 ) −1 P ≈ (Ω m 0 ) −1 P .and using the relations Also introducing the effective quantities θ(t) = θ m (t)/(p/2) and Ω 0 = Ω m 0 /(p/2) between the mechanical quantities and their electrical counterparts for a synchronous machine with p poles per phase, we obtain the more common version of the swing equation describing the dynamical relation between the rate of change of the electrical load angle and the power transmission between units: (3.7) Here M := IΩ 0 /(p/2) 2 and D := D m Ω 0 /(p/2) 2 are respectively the angular momentum of the rotor operated at the nominal grid frequency and the damping coefficient of the machine.On the right hand side of (3.7), P m denotes the net injected mechanical power (positive when the machine is operated as a generator and negative when operated as a motor), and P e denotes the electrical power injected to the grid by the synchronous machine.In the oscillator model we again assume perfect voltage support (Assumption 3.1) and lossless transmission lines (Assumption 3.2), which yield Putting everything together, we obtain the governing equations of the oscillator model of AC power grids with the parameters P i := P m i /M i , α i := Di /M i , and Proposition 3.1 (Linear stability of the oscillator model) AC power grid systems described by the oscillator model (3.8) with underlying interaction topology G(V, E) is at least neutrally stable at a fixed point θ * , if all edges are not overloaded, i.e.
Proof At a fixed point of the system θ = θ * , a small deviation of the oscillators' angles Θ := θ − θ * follows the linear dynamics where the Jacobian matrix J ∈ R 2N ×2N of the 2N -dimensional dynamical system is given by Here L is a weighted graph Laplacian as defined in (2.9) with γ ij = cos(θ * j − θ * i ), A is an N × N diagonal matrix with A ii := α i , and 0 N and I N are respectively the N × N zero matrix and identity matrix.
Let w = (w 1 , w 2 ) ∈ C 2N with w 1 , w 2 ∈ C N be an eigenvector of J corresponding to eigenvalue µ ∈ C. By definition we have J w = µw, which by writing (3.9) as a second order differential equation implies (3.11) Multiplying both sides of the equation above with the conjugate transpose, w † 1 of w 1 from the left, we obtain an expression of the eigenvalue for all (i, j) ∈ E (see Proposition 2.1 and Remark 2.2) [cf.Manik et al. (2014)].Therefore the eigenvalue always has a non-positive real part, implying that the networked dynamical system is at least neutrally stable at the fixed point.
Remark 3.2 (Neutral stability and global phase shift) Connected AC power grids described by the oscillator model (3.8) is neutrally stable at a fixed point only when the deviation Θ is a global phase shift.Because (3.11) indicates that the Jacobian eigenvalue µ = 0 only when χ 3 = 0, which implies Lw 1 = 0 and thus w 1 lies in the Laplacian eigenspace corresponding to the eigenvalue λ [0] = 0.For connected graphs there is only one eigenvector = 0, therefore the system is only neutrally stable when all nodes undergo a phase shift with the same magnitude.A global phase shift has no effects on the power flow pattern in the network since pairwise phase differences across edges remain the same.
Proposition 3.2 (LRT of the oscillator model and homogeneous nodal damping) Consider an AC power grid oscillator model with arbitrary topology (3.8) with homogeneous nodal damping α i = α ⩾ 0 for all nodes i.Then the network-wide linear responses to arbitrary external perturbations near a normal operation state θ * can be expressed explicitly in the eigenbasis of a weighted graph Laplacian: i) The network response to time-independent distributed perturbations D * is (3.12) and ii) the network response to a single sinusoidal perturbation given by D Proof Vectorizing the linear response Θ(t) of the system (3.8) to a perturbation vector D(t), we obtain the matrix equation describing the response dynamics of the oscillator model of AC power grids where α := (α 1 , • • • , α N ) denotes the vector of damping parameters and " • " denotes the Schur (element-wise) product of two vectors.Let α i = α for all i, the term α • Θ reduces to a scalar multiplication α Θ, thereby all terms involving the variable Θ in equation (3.14) can be expressed as linear combinations of Laplacian eigenvectors.Because the L here is real and symmetric so that we can write ℓ] .Using the same trick as in Remark 2.3, we obtain equations for the time-dependent coefficients c [ℓ] (t) given by (3.15) Assuming at t = 0 the AC power grid system operates normally at the fixed point θ * , we have initial conditions for the coefficients c [ℓ] (0) = 0 and ċ[ℓ] (0) = 0 for all ℓ, thereby obtain explicitly solutions for the coefficients and for the network linear responses.
For perturbations independent of time, D(t) = D * , such as constant shifts in power injection and consumption, the linear response of a power grid system is given by directly solving for the coefficients in (3.15): ℓ] and η [ℓ] := α 2 /4 − λ [ℓ] .For time-dependent perturbations D(t), such as fluctuating power injections from renewables, we obtain the network linear response based on the responses to each single frequency components at each perturbed nodes, as discussed in Remark 2.4.Similarly, we let D i (t) = δ ik εe ı(ωt+φ) and obtain the oscillator model's linear response to a sinusoidal signal at node k as (3.17) Remark 3.3 (Low dissipation regime and grid eigenfrequencies) In case the dissipation in the system is low enough such that α < 2 √ λ [ℓ] for the ℓ-th eigenvalue, in the solution of the linear response (3.13) the corresponding η [ℓ] for the same eigenmode becomes imaginary, suggesting this mode is oscillating under-damped in the power grid system with an exponentially decaying amplitude proportional to e − α 2 t , i.e., with a time constant τ = 2/α.Such intrinsic oscillation modes can also be identified by looking at the eigenvalues of the Jacobian matrix at the fixed point.Let w 1 in (3.11) be the ℓ-th eigenvector of the Laplacian L, we can see the corresponding Jacobian eigenvalue µ 4 of the power grid system.Since for a connected networked system with N nodes has N − 1 positive Laplacian eigenvalues, it also has a band of N − 1 eigenfrequencies if the dissipation is sufficiently low satisfying α < 2 √ λ [ℓ] for all N − 1 Laplacian eigenvalues.

Emerging network response patterns from LRT
In Sec. 3 we applied LRT on power grid models and obtained explicit solutions for linear network responses to perturbations at a normal operation state.The solutions are expressed in terms of the eigensystem of a weighted graph Laplacian.These Laplacians and thus their eigensystems contain information about the underlying network topology as well as the base operating state of the system, enabling us to understand and to manipulate how complex networked systems such as power grids collectively respond to external perturbation signals.
In this section we focus on the dynamic responses of AC power grids to time-varying perturbations based on the oscillator model (see Sec. 3.2) and explicate how steady-state response patterns constituted by the set of nodal response magnitudes as well as transient response patterns describing the spatio-temporal spreading of a perturbation in power grids are mathematically extracted from the solution given by the LRT.

Frequency regimes of steady-state response patterns
After a transient phase characterized by a dissipation-related time constant τ = 2/α (cf.Remark 3.3), the perturbed power grid systems reside in a second regime of network responses, where the entire network respond periodically to perturbation signals for t ≫ τ (see 3.12 and 3.13).We thus call the network responses for such large times steady-state responses.We remark that the steady-state responses here are not necessarily stationary, meaning that the nodal responses themselves can vary with time, but their characteristics, such as the amplitude and the phase of sinusoidal responses, constitute network-wide response patterns that are time-independent.
Proposition 4.1 (Steady-state response pattern for a constant perturbation) For AC power grids with arbitrary topologies (3.8), the steady-state responses to time-independent perturbations D(t) = D * , i.e. with a perturbation frequency ω = 0, near a normal operation state θ * are constituted by a homogeneous shift of grid frequency and a topology-dependent phase shift Proof By definition, the steady-state response patterns become clear by investigating the asymptotic behaviour of the responses as t → ∞.For responses to a constant perturbation vector D(t) = D * (3.12) the steady-state response reads which consists of two characteristic patterns: i) the phases of all units drift away from the normal operation state with a constant angular velocity i /α, and ii) a time-independent and unit-specific phase shift.Since v [0] = 1 √ N 1, the former pattern represents the constant homogeneous shift of grid frequency (4.1) and the latter the topology-dependent phase shift (4.2).

Remark 4.1
The global grid frequency shift (4.1) is a consequence of the imbalance between power injected into and drawn from the power grid system, which is imposed by the constant perturbation D * .The topology-dependent phase shifts (4.2) at all units suggest a network-wide redistribution of power flows on the transmission lines.The first-order approximation of the change of the load L ij := sin(θ j − θ i ) on line (i, j), δL ij = cos(θ * j − θ * i )(δθ j − δθ i ), provides an indicator for the emerging risks such as overheating for heavily-loaded lines with L ij approaching the upper limit of its safety range.Proposition 4.2 (Steady-state response pattern for a sinusoidal perturbation) For AC power grids with arbitrary topologies (3.8), the steady-state responses to a sinusoidal perturbation at node k, i.e.D(t) with D i (t) = δ ik εe ı(ωt+φ) , near a normal operation state θ * are constituted by a homogeneous phase shift and sinusoidal responses at each node with the same frequency ω and a characteristic complex amplitude for node i.

Proof
The steady-state response to a sinusoidal perturbation at a single node k is obtained by studying the asymptotic behaviour of the responses (3.13) has the form of which is composed of a homogeneous phase shift ıεe ıφ αωN and a driven oscillation at each node.Each node's angular variable θ i changes at the same frequency as the perturbation frequency ω, but with a complex amplitude (4.7) Remark 4.2 (Characterization of the steady-state response patterns to a sinusoidal perturbation) The complex nature of the amplitude suggests shifts in the amplitude and in the phase between the perturbation signal and the response signal, which are both topology-dependent and node-specific.Hence we also refer to R as the nodal response factor of node i to a sinusoidal perturbation at node k.
The homogeneous phase shift contributes to neither the change of grid frequency nor the overall power flow pattern in the network, while the absolute value of R (k) i determines the maximal deviation of the local grid frequency at node i caused by a perturbation at node k If it exceeds the safety range of normal operation, the local frequency deviation may cause damage the synchronous machine and other related grid components such as the turbine.In the rest of the subsection we focus on the steady-state response pattern constituted by the set of nodal response strengths for each node i and discuss in detail its distinctive spatial distributions in different frequency regimes.

Distance to driven unit d
ω/π = Hz < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > .The black square marks the perturbed node.Network settings are the same as Fig. 2 in Zhang et al. (2019).

Regime of grid resonances
As suggested in Remark 3.3, the dynamics of the perturbed oscillator model of AC power grids can be understood in comparison with the dynamics of a driven damped harmonic oscillator.For each intrinsic under-damped oscillation mode corresponding to a non-zero Laplacian eigenvalue λ [ℓ] > α 2 /4, the oscillatory power grid system resonates if the perturbation frequency matches the corresponding eigenfrequency ω eigen , the network responses exhibit large amplitudes since the rationalized denominator (ω 2 − λ [ℓ] ) 2 − α 2 ω 2 for the corresponding ℓ-th oscillation mode is minimized.However, the response amplitudes vary greatly for different nodes in the network due to the factor v i .If the system dissipation is sufficiently low such that there exist N − 1 under-damped oscillation modes, the corresponding N − 1 eigenfrequencies form a resonance regime where the power grid system can potentially exhibit strong distributed responses across the network (cf.Remark 3.3).
We emphasize that the spatiotemporal resonance pattern is characteristic of each perturbation frequency within the resonance regime, of each specific network topology including the prior perturbation base state, and of each location of perturbation.Therefore a power grid system's responses to real-world fluctuating perturbations containing a collection of frequency components within the resonance regime are quite irregular both temporally and spatially, which makes it a non-trivial task to evaluate the risks in perturbed power grids induced by network resonances (see Sec. 4.3 for further discussions).

Homogeneous responses: the low-frequency regime
The network response pattern for lower perturbation frequencies, i.e. the ones lower than the smallest eigenfrequency ω , can be understood by investigating the asymptotic behaviour of the nodal response strength Proposition 4.3 (Homogeneous responses at the low frequency limit) As the perturbation frequency ω → 0, the steady-state response strength at each node of an AC power grid system with an arbitrary topology (3.8) approaches a constant value, i.e.
Proof Considering (4.5) and λ [0] = 0, we derive the asymptotic behaviour of the real part and the imaginary part of the nodal response factor R As the asymptotic behaviour of the response strength i | is dominated by the imaginary part, we have Remark 4.3 (Consistency with the homogeneous grid frequency shift at ω = 0) The homogeneous response strength at each node as ω → 0 suggests a global shift of grid frequency inversely proportional to the network size N and the system dissipation parameter α.This result is quantitatively consistent with the homogeneous grid frequency shift induced by constant perturbations, as discussed in Proposition 4.1 with D * i = δ ik if node k is the perturbed one.
Remark 4.4 (Relative nodal response strength) The homogeneous nodal response at the low frequency limit (4.9) serves as a reference value for the nodal response strengths of a sinusoidally driven network.Therefore, by normalizing A (k) i with its low frequency limit 1 N α , we define the relative nodal response strength so that the nodal response strengths can be compared across networks with different sizes and dissipation values.

Localized responses: the high-frequency regime
Similarly, we investigate the network response pattern in the high-frequency regime where ω > ω eigen [N −1] by observing the asymptotic behaviour of the nodal response strengths as the perturbation frequency becomes sufficiently large and approaches infinity.where ki is a distance-and node-specific prefactor but independent on the perturbation frequency.
Proof To determine the asymptotic behaviour of the response strength i ] from (4.10) and (4.11) to a common denominator M (ω) and obtain the respective numerators N Re (ω) and N Im (ω) as polynomials of ω, (ω), and (4.16) (ω) with (4.17) The asymptotic behaviour of M (ω), N Re ki (ω) and N Im ki (ω) as ω → ∞ is dominated by the respective leading terms with the highest power of ω.The denominator scales asymptotically as For the numerators the leading term depends on a common product Q [ℓ] (ω).As shown in Appendix A, Q [ℓ] (ω) explicitly depends on λ [ℓ] in the form of Q [ℓ] (ω) can be written in terms of λ [ℓ] and other variables that are dependent on the underlying matrix L but independent of the choice of ℓ.Thus, we define where the coefficients C [j] (λ [ℓ] ) are polynomials in λ [ℓ] of degree j.Inserting the expression of Q(λ [ℓ] , ω) (4.20) to the numerators (4.16) and (4.17), we obtain where F [j] (λ [ℓ] ) and G [j] (λ [ℓ] ) are also polynomials in λ [ℓ] of degree j and can be written in terms of C [j] (λ [ℓ] ) as Considering the numerators N Re ki (ω) and N Im ki (ω) as the ki-th elements of numerator matrices N Re and N Im , we can conveniently write (4.21) in a matrix form with the coefficient matrices Φ and F [j] , G [j] are diagonal matrices with F [j] ii := F [j] (λ [i] ) and G

[j]
ii := G [j] (λ [i] ), respectively.By spelling out the polynomials m ∈ R, we can see the diagonal matrices F [j] and G [j] are polynomials of a diagonal matrix Λ with Λ ii := λ so that the coefficient matrices Φ [j] and Γ [j] can be written as g [j]  m L m , (4.28) indicating that they are polynomials of the weighted graph Laplacian matrix L of degree j, as VΛV T ≡ L. In short, the numerators N Re and N Im ki are in fact polynomials in ω, with its coefficient being also a polynomial in L, We note that as the index j increases, the powers of ω decrease and the degrees of the polynomials Φ [j] (L) and Γ [j] (L), i.e. the coefficients of ω, increase.For sufficiently large perturbation frequency ω, the leading terms in the numerators which dominate the asymptotic behaviors would be the ones with the highest degrees of ω with nonzero coefficients.We know from graph theory that (L m ) ij ̸ = 0 only for node pair (i, j) with the graph theoretic distance d(i, j) between them, i.e. the length of the shortest path from j to i on the unweighted graph defined by L, satisfying d(i, j) ⩽ m.Therefore, the first terms in N Re ki (ω) and N Im ki (ω) with the highest degrees of ω have exactly zero coefficients, i.e.Φ [j] = 0 and Γ [j] = 0, for all j < d(k, i).Consequently, the leading term in the numerators are ki ω 4N −3−2d (4.30) with d := d(k, i).Together with the asymptotic behavior of the denominator (4.19), we obtain ki ω −2d−1 .(4.31) Remark 4.5 (Localized response patterns in the high-frequency limit) Proposition 4.4 implies that the response amplitude of the grid frequency A (k) i decays exponentially with the graph-theoretic distance d between the perturbed node and the responding node.The response amplitude also decays as a power law in ω for fixed d.In the high-frequency limit, a network's response to the perturbation is restricted to the perturbed node, i.e.
with δ ki being the Kronecker delta function.
namics.Each unit i has n state variables (x ) and is governed by dynamics where f : R N → R and g : R → R are functions respectively representing the intrinsic and the coupling dynamics of units and allowing for a stable fixed point of the system.If unit k is sinusoidally driven with frequency ω → ∞, we conjecture that the amplitude Ã(k) i,m of the sinusoidal response in state variable D m t x i of unit i is given by (4.33) where ki | is a distance-and node-specific prefactor but independent on the perturbation frequency.
The response of a sinusoidally driven damped harmonic oscillator with Ã(k) i,m ω→∞ ∼ ω −2 can be seen as a special case of (4.33) with n = 2, d = 0 and m = 0.For networks of Kuramoto phase oscillators with n = 1, (4.33) is proven to be valid [see the Supplementary of Zhang et al. (2019)].For the oscillator model of power grid with n = 2 and m = 1, (4.33) reduces to Proposition 4.4.
Remark 4.7 (Generalizability of the steady-state response patterns in three frequency regimes) In the above discussions of the steady responses patterns in three frequency regimes in Sec.4.1, we do not make any assumptions on the network topology.Therefore, our results on the characteristics of the homogeneous, the resonance and the localized response patterns in three regimes hold for arbitrary network topologies.Nevertheless, the evaluation of the parameters and the prefactors, such as Φ [d] ki in (4.14), is network-and topology-dependent by definition.

Topological factor of transient spreading dynamics
In the following we focus on the transient response of AC power grids to external perturbations, which is referred to as the network responses close to the time of perturbation and thus describes the spatiotemporal pattern in the perturbation spreading process across the network.Particularly, we demonstrate how to extract the role of the network topology in the spreading pattern based on the LRT of the oscillator model.
To investigate the transient response close to the onset of perturbation at t = 0, we Taylorexpand the linear response (3.13) at node i to a sinusoidal perturbation at node k in powers of t as around t = 0, which is characterized by the time derivatives of the linear response at t = 0.
Here D n t := d n dt n is Euler's notation for differential operator.The n-th order time derivative of the linear response at t = 0 is (4.35)The transient response of a power grid network can thus be estimated by the first non-zero term in the power series of t (4.34).Interestingly, the resulting series does not start with low powers of t such as t 0 or t 1 as typical for common Taylor expansions.Instead, it typically starts with large powers of t as the following proposition illustrates.Proof To find out the leading term in the Taylor series of the linear response (4.34), we study the summand of the ℓ-th eigenmode in the derivative of the linear response i εe ıφ F n (λ [ℓ] ) for convenience5.In the summand, the function F n (λ [ℓ] ) appears to be a division of two polynomials of λ [ℓ] .In fact, it can be shown that the leading term of F n (λ [ℓ] ) (denoted as LT F n (λ [ℓ] ) ), i.e. the term with the highest order of λ [ℓ] is if n is even. (4.37) A derivation of the result (4.37) is given in Appendix B. We note that, for n = 0 and n = 1, F n (λ [ℓ] ) = 0, which is a consequence of the choice of initial condition: the linear response and its first derivative are supposed to be zero at t = 0, as the responses and the frequency response Θ(k) i are zero at the onset of perturbation.For n ⩾ 2, (4.37) indicates a monotonic relation between the degree of F n (λ [ℓ] ) as a polynomial of λ [ℓ] and the order of derivative n: As we have shown in Sec.4.1.3,sums of the form of ), where P j (λ [ℓ] ) represents a general polynomial of λ [ℓ] of degree j, can be seen as [P j (L)] ki , the ki-th element of the matrix P j (L), a polynomial of L with degree j.Applying this result to the derivatives of the linear response (4.35), we find that D n t Θ (k) i (0) can be considered as the ki-th element of matrix F n (L), a polynomial of L. The leading term of D n t Θ (k) i (0) is thus given by (4.38) which contains (L m ) ki with m = n−3 2 if n is odd and m = n−2 2 if n is even.We notice that for a given node pair (k, i) at distance d, we have (L m ) ki = 0 for all m < d because no path of length m < d can connect nodes k and i. Therefore all terms in the Tayler series (4.34) with the leading term's degree lower than d vanish.The first non-zero term in the series thus equals the leading term of the (2d + 2)-th derivative of the linear response, because all other terms contain (L m ) ki with m < d and thus vanish.Taken together, the transient linear response near t = 0 can be approximated as Remark 4.8 (Topological factor in perturbation spreading) The leading-term approximation of the linear response (Proposition 4.5) provides a way to disentangle the impact of various factors on the dynamical spreading process in power grid networks.Specifically, the impact of the specific network topology, including the interaction structure between units and the system's base state at t = 0, is reflected in the factor (L d ) ki in the leading-term approximation.It satisfies (4.41) suggesting that it can be interpreted as the product of the edge weights along a shortest path P d k→i between node k and i, summed over all shortest paths.This insight provides guidelines for manipulating the perturbation spreading dynamics in power grid networks through changing the underlying topology.Numerical evidences show that the topological factor that revealed by the leading-term approximation also enables a master function approach to accurately predict threshold-crossing arrival times in power grid networks [Zhang, Witthaut and Timme (2020)].
Remark 4.9 (Scaling behaviours in transient responses) The leading-term approximation of the transient response (4.36) reveals two scaling behaviours as t → 0 (cf.Fig. 3).First, the transient response grows algebraically in time with a distance-dependent exponent: ! is a time independent prefactor but depends on signal magnitude, topology, base operating state and inter-node distance d.Second, the transient response decays nearly exponentially with distance d, since the factor t 2d+2 dominates the asymptotic behaviour of the response at large but finite distances as t → 0. < l a t e x i t s h a 1 _ b a s e 6 4 = " as t → 0, with time independent constant C d that depends on signal magnitude, topology, base operating state and inter-node distance d, see (4.36).Thus, responses (b) algebraically increase with time t at any given unit and (c) at any given time, they decay nearly exponentially with shortest-path distance d = d(k, i) between the perturbed unit k and the observed unit i.The grey dotted lines in (b) indicate the leading-term approximations.Network settings: (N, N g , P g , P c , K g , K c , α) = (6, 3, 1 s −2 , −1 s −2 , 10 s −2 , 10 s −2 , 1 s −1 ).For the perturbation signal (ε, ω/2π, φ) = (1, 1 Hz, 0 rad).
Remark 4.10 (Generalizability of the transient spreading patterns) The above discussions on the spatio-temporal pattern of transient spreading do not involve any assumptions on the underlying network topology.Thus the form of the leading term approximation of the transient network response (Proposition 4.5) does not depend on the specific choice of network topology.We underline that the evaluation of the topological factor in perturbation spreading (Remark 4.8) is indeed topology-and node-specific, as it captures the local interaction structure and consists of all shortest paths between the perturbed node and the responding node.

Nodal vulnerability to unpredictable fluctuations
In Sec.4.1 and Sec.4.2 we show how the distributed response patterns of power grid systems, in the steady-state (t → ∞) and in the transient stage (t → 0), are analytically extracted from the LRT of the oscillator model.The response patterns are numerically proven to be highly accurate [Zhang et al. (2019), Zhang, Witthaut and Timme (2020)], but the results are valid only for given perturbation signals, hence deterministic.Meanwhile, real-world power grid systems are perturbed by power fluctuations whose exact time series are hardly predictable.In this section we discuss how LRT helps to estimate network responses to random perturbations.
As discussed in Sec.4.1, power grid systems respond resonantly to perturbations with frequencies falling in the band of network's eigenfrequencies eigen ], exhibiting the most irregular network-wide spatiotemporal patterns, compared to the almost homogeneous pattern for lower frequencies and the localized pattern for higher frequencies.Moreover, the resonance response pattern varies drastically for different perturbation frequencies and for different locations of perturbation.Therefore, estimating the nodal responses for random perturbation signals involving frequency components within I res is a task not only of practical significance regarding the operational safety of power grid systems, but also with a high theoretical complexity.eigen ], the nodal Dynamic Vulnerability Index (DVI) is defined as Proof To analyse the steady network response, i.e. the response in a time window T → ∞, to a perturbation signal composed of a range of frequencies, we make use of the steady-state response of the oscillator model to a sinusoidal signal (4.6).The steady-state response Θi (t) of the frequency at node i to a perturbation signal εe ı(ωt+φ) at node k is given by ωt+φ) , (4.44) which can be seen as a driven oscillation with a complex amplitude ıεωR As a consequence of the linear nature of LRT, the nodal frequency response to a temporally fluctuating perturbation signal containing a spectrum of frequency components is obtained by summing the response Θ(k) i (ε, ω, t) (4.44) over all frequency components (see also Remark 2.4).Particularly, for perturbation signals which are characterized by a specific power spectral density (PSD) S(w), the strength of its frequency component εe ı(ωt+φ) can be expressed in terms of the frequency as ε(ω) ∝ S(ω) 1 2 while no assumptions is made on the choice of its phase φ.For instance, in modern power grids integrated with renewable energies, the power fluctuations from wind and solar energy are characterized by a power law PSD with the Kolmogorov exponent −5/3 [Anvari et al. (2016)].For such resonant perturbations, the nodal frequency response reads ωt+φ) dω. (4.45) Here c is a scaling factor between the power distribution over frequencies of a specific signal and the (normalized) PSD.Independent of the specific realization of the fluctuation signal, the largest possible magnitude of the frequency response (4.45) is reached only when the involved frequencies are finite in number and non-resonant to each other6 so that all oscillating responses for each frequency component ω (4.44) would eventually align.The largest possible magnitude of frequency response time series is given by the sum over ω of the magnitude of each frequency response |ıcS(ω) For finite time series with M data points and time step ∆T , its frequency components given by discrete Fourier transform are apparently resonant to each other.However, for time series with a fixed sampling rate, the longer the observation time window T , the smaller the frequency interval ∆ω = 2π/T and thus the more frequency components exist within the given interval of interest I res .As T approaches infinity, the order of the resonant frequencies {ω n } becomes sufficiently high so that the alignment of the phases can be attained at a finite rate (?).Therefore, as T → ∞, the all-time-high frequency response approaches the sum of the frequency-specific response magnitudes, which is proportional to the DVI given in (4.42): If one only considers the relative ranking of the all-time-high frequency response of nodes in a given network, but not their absolute values, the overall scaling factor c in (4.46) does not play a role in the ranking and we finally arrive at lim Remark 4.11 (Generalization of DVI and convergence time) The DVI defined in (4.42) provides a measure to estimate the relative nodal risk from the power grid network resonances induced by a unpredictable perturbation signal and thus helps to identify the most vulnerable nodes in networks with arbitrary topologies exhibiting particularly strong resonant responses.The integration interval of the frequencies in DVI is chosen to be the band of eigenfrequencies I res for a specific network so that its irregular resonance patterns are covered, however it is not a must.In principle, any frequency range can be chosen for DVI to estimate the relative strength of the all-time-high responses in the specific frequency range.However, one should note that the timescale for the ranking of all-time-high responses to converge to the ranking given by DVI (in 4.46 and 4.47) depends on the chosen frequency range.For instance, the convergence time would be longer if lower frequencies are included in the integration interval of DVI.

Role of LRT in uncovering response patterns
In the previous section we elaborated how network-wide dynamic response patterns of power grid systems can be extracted from the explicit solution of nodal responses given by the LRT of the oscillator model (see Sec. 3).In this section we summarize and compare the role of LRT in revealing different categories of the dynamic response patterns, such as the patterns emerging on different timescales, and in responses to perturbations with different levels of randomness and magnitude.

Transient vs. steady-state responses
As discussed in Sec.4.1 and Sec.4.2, power grid transmission networks exhibit distinctive spatiotemporal response patterns on different timescales.The transient spreading pattern (see Sec. 4.2) of an external perturbation signal in a normally operating power grid network appears close to the time of impact t = 0.It is characterized by a set of points in time, at which the impact arrives at individual units in the network.To a large extent, the topological dependence of the arrival times can be captured by a topological factor which arises from the leading-term approximation of the linear response.The steady-state response patterns, in contrast, emerge as t → ∞ (see Sec. 4.1), where the nodal responses to a sinusoidal perturbation converge to sinusoidal oscillations as well, but with various amplitudes.Consequently the set of the nodal response amplitudes, a time-invariant but frequency-dependent feature of the oscillating responses, constitute the steady-state response patterns characterizing three frequency regimes.
The time scale separating transient from steady-state regimes is not a universal constant but intricately depends on several factors including network size N and network topology, damping constant α, and the specific node location we are interested in within the network.Both transient and steady-state response patterns have been revealed and characterized through asymptotic analyses of the explicit solution of the linear nodal responses (3.13).The solution depends explicitly on time while the dependence on the network topology is implicit, embedded in the eigenvalues and eigenvectors of the weighted graph Laplacian matrix L. Through asymptotic analyses, either with respect to time t or the perturbation frequency ω, the lengthy solution of the linear nodal response is reduced to one term per eigenmode that dominates the asymptotic behaviour as the variable t or ω approaches its corresponding limit (see Proposition 4.4 and Proposition 4.5).As the contribution of each eigenmode contains the "overlap factor" v k of the perturbed node k and the responding node i, the powers of the Laplacian eigenvalue (λ [ℓ] ) m that is involved in the dominating term in each eigenmode ℓ translates to elements of the power of the Laplacian matrix (L m ) ki through the summation over all N eigenmodes ℓ ∈ {0, • • • , N − 1}.Furthermore, with the help of the result from graph theory that (L m ) ki ≡ 0 if m ∈ N 0 is smaller than d(k, i), the shortest path distance between node k and i, the dependence of the nodal linear responses on graph-theoretic distance emerges.In this way one obtains the asymptotic spatiotemporal response patterns that depends explicitly on distance.
A major difference between the patterns we uncover in steady-state responses and in transient responses is, the asymptotic behaviors of the steady-state response patterns are exact, in the sense that the higher order terms are negligible at the observed limits of the perturbation frequency ω (see Sec. 4.1).Meanwhile, the contribution of the higher order terms are not negligible in the patterns in transient responses: the leading-term approximation of the response exhibits a diverging error as the perturbation spreads further and the arrival time grows larger.Numerical simulations show that the higher order terms accounts for about 10% of the actual arrival time of perturbations [Zhang, Witthaut and Timme (2020)], which is significant in predicting the perturbation spreading behaviour in real-world power grid systems.However, by means of numerical techniques, we can still use the topological factor proposed in Remark 4.8 to estimate the contribution of higher order terms O(2d + 3) in the Taylor series (4.36) in a specific network ensemble and give accurate predictions for the actual arrival times of the impact of a perturbation [Zhang, Witthaut and Timme (2020)].Related recent works Wolter et al. (2018), Schröder et al. (2019) studied transient propagation of perturbations in networked systems consisting of one-dimensional dynamical units.One main finding is a similar scaling of the unit's state variables x i (t) (or their deviations from a base state) with time t as x i (t) ∼ t d where d is the shortest-path distance between perturbed node and the node i the response in measured at.

Responses to deterministic perturbations vs. responses to unpredictable perturbations
The power grid response patterns we discuss in this article can be classified into two categories: the ones emerging in the deterministic responses to a given perturbation signal, such as the transient responses (Sec.4.2) and the steady-state responses (Sec.4.1) to a given signal, and the ones that estimate the cumulative impact of an unpredictable signal on a power grid network, such as the DVI measuring the nodal risk of network resonances (Sec.4.3).Both categories of power grid response patterns are discovered based on the explicit solution of the linear nodal responses (3.13) that derived from the LRT.
Looking more closely, one finds that the estimated patterns in the cumulative responses can be seen as a result built upon the deterministic steady-state linear responses with one more dimension, i.e. the specific properties of the perturbation signals.At its core, the steady-state network frequency responses Θ (4.44) to a single-frequency sinusoidal oscillation D k (t) at a given node k is a mapping f G,θ * : R → R N with f G,θ * depending on the underlying network topology G and network's base state θ * prior to perturbations.The fluctuating nature of the perturbation signal adds another dimension to the responses: the amplitude of D k is no longer a given constant ε, but becomes further dependent of the frequency ω through the PSD S(ω) of the signal, i.e. ε(ω) = S(ω) 1 2 .In this way, the unknown or unpredictable temporally detailed features of the perturbations are integrated into the LRT framework for networked dynamical systems, which extends standard LRT statements and enables us to estimate features of network responses beyond the deterministic realm of systems driven with known signals.
We emphasize that, due to the intrinsic irregularity of fluctuating perturbation signals, the errors of the estimates for the associated responses appear to be significantly higher than the ones for the deterministic responses [compare Zhang et al. (2019) and Zhang, Ma and Timme (2020)].In a finite signal time series characterized by a PSD function, the contained frequencies are finite in number in the considered frequency interval (such as the resonance regime I res ), and apparently resonant to each other, which leads to a deviation in the ranking of the all-time-high nodal responses to the ranking given by indices computed a priori (such as the DVI discussed in Sec.4.3).Additionally, realistic perturbation signals do not follow exactly the characteristic PSD, such that randomness also exists in the amplitudes of the frequency components.Nevertheless, compared to the deterministic patterns which gives only a posteriori information of network responses, estimates such as the DVI given by the extended LRT may provide a useful guiding tool for risk assessments in real-world power grid systems.

Small responses vs. large responses
As the name suggests, LRT provides the linear approximation of a system's response to a perturbation close to a fixed point of a networked dynamical system.Therefore the solution given by LRT intrinsically deviates from the actual system responses due to the neglected higher order terms in the system's collective nonlinear dynamics.As the system being driven further and further away from the fixed point, the responses typically increase and so do the estimation error of the LRT.However the range of validity of the LRT, as well as how its error grows with the perturbation, is usually nontrivial and system-dependent.
For the oscillator model of power grid networks, the error of the solution given by LRT (3.12 and 3.13) follows the same trend and grows with an increasingly stronger perturbation signal.However, numerical evidence shows that the error increases mildly with the magnitude of the perturbation until it blows up close to a bifurcation point [Zhang et al. (2019)].In the linearized dynamics of the oscillator model at a fixed point θ * (3.14), the deviation of the nonlinear coupling terms sin(θ j − θ i ) for all edges (i, j) ∈ E to their values sin(θ * j − θ * i ) at the fixed point are represented by the first-order approximations cos(θ * j −θ * i )(Θ j −Θ i ), vectorized as the term −LΘ in (3.14).For power grid systems working at a stable operation state without any transmission line overloaded (|θ * j − θ * i | ⩽ π 2 for all edges (i, j) ∈ E, see Proposition 2.1), the linear approximation breaks down only when the system is driven far enough from the fixed point θ * and goes close to the point where one of the lines (i, j) is fully loaded, i.e. sin(θ j − θ i ) = 1.In this regime, the linear approximation diverges from the actual state of the system and the error grows explosively.As elaborated in an article by Manik et al. (2014), when one of the transmission line is fully loaded, the power grid system reaches a bifurcation point where the initial stable fixed point is lost and the oscillating units in the system becomes desynchronized.Therefore, the LRT of the oscillator model of the power grid systems, together with all of the derived response patterns, are generally valid as long as none of the lines become overloaded and the entire system becomes unstable [see Zhang et al. (2019) for quantitative results of the LRT errors].

Conclusions and discussions
In this work we systematically discuss how linear response theory (LRT) may shed light on the spatio-temporal response patterns emerging in networked dynamical systems under timedependent perturbations.We exemplify a full analysis for model dynamics of power grid systems which are inevitably exposed to fluctuating power injections from renewable energy sources.Beyond previous works, we integrate and present all details required for a full mathematical analysis, specifically demonstrate how to evaluate the generally intricate, multiple-sum expressions determining spatio-temporal response patterns in a useful way and highlight how different results interconnect, for instance between transient and long-term dynamics or between different types of perturbations and across topologies.We introduce the main ideas of LRT and its general requirements for applicability on system settings (cf.Sec. 2).We explicate various aspects of application to models of power grid systems, such as i) the solution of linear responses of the stationary DC power flow model and of the dynamic oscillator model of AC power grids (cf.Sec.3), and ii) how it helps to identify dynamic patterns in network-wide responses which provide theoretical guidelines for power grid design, control and risk assessments.
Although LRT has been widely used as a powerful tool in analyzing various response dynamics of many complex networked dynamical systems, the works presented in this article provide a fresh methodological angle to approach the problem.For power grid systems, LRT has been used to estimate e.g.quadratic performance measures for the network's overall excursion away from synchrony [Tyloo et al. (2018), Tyloo et al. (2019), Plietzsch et al. (2019), Coletta et al. (2018)] and the variance of the frequency response increment distribution [Haehne et al. (2019)].For connectome dynamics in brain, LRT has also been used to analytically estimate the covariance of Gaussian linear model of the stochastically perturbed system [Tononi et al. (1994)], which links the structural and the functional connectivities between brain regions [Zamora-López (2010), Wang et al. (2019)].
The work presented in this article approaches the dynamic network responses in a way different from the above-mentioned works: instead of quantifying the stochastic features of the overall or distance-specific responses directly based on the linear responses, we start from explicating the deterministic solution of network-wide responses to a single-frequency signal and using methods from graph theory and asymptotic analysis to extract spatiotemporal response patterns.These may be interpreted in a physically intuitive way.Specifically, the three frequency regimes of steady-state response patterns (Sec.4.1) and the master curve of transient perturbation spreading (Sec.4.2) are entirely deterministic.Especially, the former work on fluctuation-induced network resonances can be seen as a direct generalization of the classical resonance phenomenon of a single driven damped harmonic oscillator to oscillators interacting on networks.The emerging pattern constituted by the estimations of the all-time-high nodal response magnitudes to a irregularly varying perturbation (Sec.4.3), is also a straightforward result derived from the spatial patterns in the network-wide responses to a deterministic (periodic) signal.

Challenges and future work
Future work regarding the response theory for networked dynamical systems may follow several directions.
First, exact analytical solutions, and even many asymptotic results of the linear responses of general networked dynamical systems, as well as the response patterns emerging from these solutions, remain unknown to date.Applying the LRT presented in this work on a networked dynamical system, we employed several conditions on the system's dynamics (see Sec. 2) to explicate a full analysis without too many notational and other complications.The conditions include homogeneous nodal dynamics, a diffusive coupling term g ij (x j − x i ), and the evenness of the coupling function's sensitivity dgij d(xj −xi) to small changes in the difference of nodal states x j − x i , such that a symmetric weighted Laplacian matrix arises in the linearized response dynamics of the system at the fixed (operating) point (2.10).The presence of a symmetric Laplacian matrix in the linearized dynamics ensures the option to express the linear responses in the Laplacian eigenbasis, which plays a critical role in linking the response at a specific node to the graph-theoretic distance to the perturbation.Thereby it is critical also in uncovering the topological structure of the dynamical response patterns across the network.However, for many networked dynamical systems, such as the third-order model of power grid dynamics including the voltage dynamics [Machowski et al. (2008)], such preconditions are not fulfilled.One way to overcome this theoretical barrier and to extend LRT to such networks of dynamical systems is to transform the system's state variables to another coordinate system where the Jacobian matrix J in (2.3) is diagonal or almost diagonal (such as in the Jordan normal forms of J ).In this way explicit solutions for linear responses can be obtained in the new coordinate system where dynamic response patterns can be identified in similar ways presented in this work.
Second, one could use LRT to develop strategies to control the impact of fluctuations on networked dynamical systems such as power grids.So far, we gained insights into the spatiotemporal structure of the responses across networks and developed indices to estimate the nodal risks against external fluctuations.The next step towards more reliable and more robust power grid systems would be to utilize the obtained understanding to develop countermeasures against the risks, e.g. to suppress the potentially dangerous responses such as network resonances and to slow down the spreading of the impact of a sudden drop of injected power.A potential way to achieve such tasks could be to manipulate discovered response patterns by changing the interaction structure of the power grid.
Third, research on how different classes of network topologies potentially impact response patterns through their specific characteristics of their eigensystems.Progress in this direction seems hard, because one would need to be able to characterize, e.g., eigenvectors of graph ensembles such that they directly help to extract useful information from complex expressions like (4.5) or even (3.13)specifically for that ensemble.
Fourth, the LRT per se could be extended by considering also the higher-order approximations of the system's responses, cf.[Thümler et al. (2022)].In the current work, we demonstrated that many features of collective response patterns of networked dynamical systems with nonlinear couplings, such as the oscillator model of power grids, are dominantly captured by the first-order (i.e.linear) approximation of the system's responses, yet nonlinear effects may also play a role for other systems with certain forms of the intrinsic nodal dynamics or the interaction dynamics between nodes, specifically if we ask for the loss of solutions near operating states [Thümler et al. (2022)].Therefore it would be desirable if the contributions of the higher-order approximations of the system's responses can be estimated.An open question also here is how such nonlinear effects depend on the interaction topology of the network.
We conjecture that general network dynamical systems, also beyond power grids, similarly respond in characteristic ways to external input signals, making the systems non-equilibrium and often non-stationary, and to be described by non-autonomous deterministic or stochastic evolution equations.Several of the analysis steps presented above hint that the key methodological tools are either readily transferable to more general systems' settings or may be adapted to such settings.Candidate classes of systems include networks of multi-dimensional units, with discrete or hybrid dynamics, with delayed interactions or with spatially or temporally correlated stochastic inputs.Application areas may range from gene regulatory networks and metabolic circuits in cell biology to the controlled self-organized dynamics of engineered systems with feedback, from complex mechatronic systems to swarms of autonomous aerial vehicles.1,0 = Y 1,0 − (α 2 − 2λ [ℓ] ), which is a polynomial of λ [ℓ] with degree 1 since Y 1,0 is a constant independent of λ [ℓ] .n,0 (λ [ℓ] )] = n.
Step 3: Finally, we show that the sum of coefficients involving both r 2 and r 3 over s 2 ∈ S ℓ b and s 3 ∈ S ℓ \{s2} c , i.e. all possible partitions P ℓ (a, b, c), is a polynomial of λ [ℓ] with degree j = b + 2c.That is, equation (A 7).

Proposition 4. 4 (
Localized response patterns in the high-frequency regime) As the perturbation frequency ω → ∞, the steady-state nodal response strength A (k) i in an AC power grid system with an arbitrary topology (3.8) decays as a power-law of ω with an exponent depending on the graph-theoretic distance d between node k and i, i.e.
e g 7 + o F + 9 t Z 7 b 3 p B b + 9 G u r o y y 9 l C C 6 N 3 8 B s l n n K H < / l a t e x i t >
Proposition 4.5 (Leading-term approximation of transient response) The transient response at node i in an AC power grid network (3.8) to a sinusoidal perturbation of frequency ω at node k with an onset at t = 0 is approximated by the (2d + 2)-nd term in the Taylor expansion of the linear response (3.13) around t = 0Here d := d(k, i) denotes the graph-theoretic distance between node k and node i.

Figure 3 .
Figure 3. Transient Network Response Dynamics exhibits algebraic growth with time and exponential decay with shortest-path distance.(a) Basic network of N = 6 units illustrates (b,c) transient algebraic responses [color-coded as units in (a)] to a sinusoidal perturbation at node 1 that increase like Θ

Definition 4. 1 (
Dynamic vulnerability index (DVI) for random network resonances) In an AC power grid system (3.8)perturbed by a random fluctuation at node k, characterized by a power spectral density S(ω) with frequency components ω ∈ I res = [ω DVI estimates ranking of the nodal all-time-high steady frequency responses) Let an AC power grid system (3.8)be perturbed by a time-dependent fluctuation at node k.A random signal time series is characterized by a power spectral density S(ω), i.e., the strength ε of its frequency component εe ı(ωt+φ) is frequency dependent and follows ε(ω) ∝ S(ω) 1 2 , and the corresponding phase φ is randomly drawn from the uniform distribution on [0, 2π), independently for each realization of the fluctuation time series.Suppose the fluctuation signal is composed of frequency components with ω ∈ I res = [ω the ranking σ ATH of the nodal all-time-high steady frequency response magnitude max t∈[0,T ] | Θ(k) i (t)| in an observation window T approaches the ranking σ DVI of the DVI defined in Definition 4.1 as the observation time window goes to infinity, i.e. lim T →∞ σ ATH (i) = σ DVI (i) for all i ∈ {1, • • • , N }.(4.43) defined in (4.5) characterizing the response strength at individual nodes in a network.The complex amplitude gives rise to an amplitude shift |ıεωR (k) i | and a phase shift arg(ıεωR (k) i ), which both are specific to the perturbation strength ε and the perturbation frequency ω.