eurm10 msam10 remarkRemark[section] assumptionAssumption[section]
definitionDefinition[section] conjectureConjecture[section] Fluctuation
Response Patterns of Network Dynamics – an Introduction–References

# Fluctuation Response Patterns of Network Dynamics – an Introduction

XiaozhuZHANG$\,{}^{1,2,3}$ MarcTIMME$\,{}^{3,4}$ 1MOE Key Laboratory of
Advanced Micro-Structured Materials and School of Physics Science and
Engineering, Tongji University, Shanghai 200092, P. R. China
2Frontiers Science Center for Intelligent Autonomous Systems, Tongji
University, Shanghai 200092, P. R. China
3Chair for Network Dynamics, Center for Advancing Electronics Dresden (cfaed)
and Institute for Theoretical Physics, Technical University of Dresden, 01062
Dresden, Germany
4Lakeside Labs, Lakeside B04b, 9020 Klagenfurt, Austria
email: xiaozhu.zhang<EMAIL_ADDRESS>
(2021)

###### Abstract

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
for fluctuation-driven networked dynamical systems. The derivations presented
not only provide approximate analytical descriptions of the dynamical
responses of networks, but more importantly, allows to extract key qualitative
features about spatio-temporally distributed response patterns. Specifically,
we provide a general formulation of a linear response theory for perturbed
networked dynamical systems, explicate how dynamic network response patterns
arise from the solution of the linearized response dynamics, and emphasize the
role of linear response theory 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.

###### keywords:

networked dynamical systems, linear response theory, fluctuations, dynamic
responses, pattern formation, perturbation spreading, oscillatory networks,
power grid dynamics

###### 2020 Mathematics Subject Classification:

P

††volume: 000

rimary: 34C15, 94C15; Secondary: 93A14, 34E10

##  1 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.

##  2 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.

###  2.1 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
$\bm{x}=(x_{1},\cdots,x_{N})\in\mathbb{R}^{N}$. The autonomous dynamics of the
$N$-dimensional dynamical system is governed by

$\dot{\bm{x}}=\bm{f}(\bm{x}),$ (1)

where $\bm{f}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ is a function that in
general depends on the states $\bm{x}$ of all units and does not explicitly
depend on time. Let us consider a system that exhibits a fixed point at
$\bm{x}=\bm{x}^{*}$ where $\bm{f}(\bm{x}^{*})=\bm{0}$, and is influenced by an
external dynamic driving vector $\bm{D}(t)\in\mathbb{R}^{N}$ at the fixed
point $\bm{x}^{*}$. We investigate how the system dynamically responds to
$\bm{D}(t)$, i.e. how the system’s deviation $\bm{X}(t):=\bm{x}(t)-\bm{x}^{*}$
from the fixed point evolves in time. In general, the dynamics of the system’s
response $\bm{X}$ follows

$\dot{\bm{X}}=\bm{f}(\bm{x}^{*}+\bm{X})+\bm{D}(t),$ (2)

where both functions $\bm{f}$ and $\bm{D}$ can be highly nonlinear so that an
exact analytical solution of the system’s response $\bm{X}$ typically does not
exist or is unknown.

Important information about the response dynamics (2) is given by the
stability operator at the fixed point obtained from the linearization of the
function $\bm{f}$ at $\bm{x}=\bm{x}^{*}$, i.e. the Jacobian matrix
$\mathcal{J}$ with $\mathcal{J}_{ij}:=\frac{\partial f_{i}}{\partial
x_{j}}|_{\bm{x}=\bm{x}^{*}}$. The system’s response dynamics thus
approximately follows the linearized differential equation

$\dot{\bm{X}}=\mathcal{J}\bm{X}+\bm{D}(t).$ (3)

The signs of the real parts of the Jacobian eigenvalues $\bm{w}^{[\ell]}$
(with $\ell$ being the index) indicate whether the fixed point $\bm{x}^{*}$ is
linearly stable ($\mbox{Re}[\bm{w}^{[\ell]}]<0$), unstable
($\mbox{Re}[\bm{w}^{[\ell]}]>0$) or neutrally stable
($\mbox{Re}[\bm{w}^{[\ell]}]=0$) in the eigendirection/eigenspace
corresponding to the eigenvalue $\bm{w}^{[\ell]}$. Below, we typically focus
on the dynamics near stable fixed points, where all
$\mbox{Re}[\bm{w}^{[\ell]}]<0$ (or marginally stable ones where one
$\mbox{Re}[\bm{w}^{[\ell]}]=0$ if $\mathcal{J}$ is a Laplacian operator) to
justify the assumption that the solution $\bm{X}(t)$ of the linearized
equation (3) reasonably well approximates the full solution of (2).

Let us now focus on systems with pairwise interactions between the $N$ units
such that the function $f_{i}(\bm{x})$ controling the time evolution of unit
$i$ can be written as

$f_{i}(\bm{x})=h_{i}(x_{i})+\sum_{j=1;\ j\neq
i}^{N}K_{ij}g_{ij}(x_{i},x_{j}),$ (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\neq j$. Here
$K_{ij}\in\mathbb{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})\not\equiv 0$.

definition

###### Remark 2.1

Pairwise interactions (4) induce an interaction structure of the dynamical
system (1) that can be represented by a graph $G(V,E)$ where the set of
vertices $V$ consists of $N$ variables $x_{1},\cdots,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)\in V^{2}|\ (i\neq j)\land(K_{ij}\neq 0)\\}$.

The linearization of the pairwise interaction (4) yields the Jacobian matrix
of the response dynamics (2) of the driven networked system J_ij &=
∂∂xj(h_i(x_i)+∑_k≠i K_ikg_ik(x_i,x_k))—_x=x^*
= {dhidxi—_x=x^*+∑k≠iKik∂gik∂xi—x=x*for j=iK_ij∂gij∂xj—_x=x^*for j≠i. It is
clear that the topology of the interaction network $G$ explicitly enters the
Jacobian matrix through the matrix of the coupling strength
$K\in\mathbb{R}^{N\times 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 (1) near a fixed point (3) and the arising Jacobian
matrix that explicitly depends on the system’s topology (2.1) due to pairwise
interactions (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
(3), $X(t)$, can be characterized by evaluating the specific spectral
properties of the Jacobian matrix $\mathcal{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.

###  2.2 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 (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 $\tilde{g}_{ij}$ 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})=\tilde{g}_{ij}(x_{j}-x_{i})$ in (2.1). For notational
simplicity, we again denote the functions $\tilde{g}_{ij}$ just by $g_{ij}$.

######  Proposition 2.2 (Stability of diffusively coupled networks: a special
case)

A networked dynamical system with evolution function (4) and interaction
network $G(V,E)$ is at least neutrally stable at its fixed point
$\bm{x}=\bm{x}^{*}$ if a) the intrinsic nodal dynamics $h_{i}$ satisfies
$\frac{\mathrm{d}h_{i}}{\mathrm{d}x_{i}}|_{\bm{x}=\bm{x}^{*}}\leq 0$ for all
nodes $i\in V$, b) the coupling function $g_{ij}$ is diffusive and satisfies
$\frac{\mathrm{d}g_{ij}}{\mathrm{d}(x_{j}-x_{i})}|_{\bm{x}=\bm{x}^{*}}\geq 0$
for all node pairs $(i,j)\in E$.

###### Proof 2.3

The diffusive form of the coupling function $g_{ij}(x_{j}-x_{i})$ yields a
useful relation $\frac{\partial g_{ij}}{\partial x_{j}}=-\frac{\partial
g_{ij}}{\partial x_{i}}=\frac{\mathrm{d}g_{ij}}{\mathrm{d}(x_{j}-x_{i})}$.
With this particular symmetry, the Jacobian matrix (2.1) of the system at the
fixed point $\bm{x}=\bm{x}^{*}$ takes the following form

$\mathcal{J}_{ij}=\left\\{\begin{array}[]{ll}-\beta_{i}-\displaystyle\sum_{k\neq
i}K_{ik}\gamma_{ik}&\mathrm{for}\quad i=j\\\
K_{ij}\gamma_{ij}&\mathrm{for}\quad i\neq j,\end{array}\right.$ (5)

where
$\beta_{i}:=-\left.\frac{\mathrm{d}h_{i}}{\mathrm{d}x_{i}}\right|_{\bm{x}=\bm{x}^{*}}$
and
$\gamma_{ij}:=\left.\frac{\mathrm{d}g_{ij}}{\mathrm{d}(x_{j}-x_{i})}\right|_{\bm{x}=\bm{x}^{*}}$.
Given that $\beta_{i}\geq 0$, $\gamma_{ij}\geq 0$ and $K_{ij}\geq 0$ by
definition, the Jacobian $\mathcal{J}$ is diagonally dominant:

$\left|\mathcal{J}_{ii}\right|=\left|-\beta_{i}-\displaystyle\sum_{k\neq
i}K_{ik}\gamma_{ik}\right|=\beta_{i}+\displaystyle\sum_{k\neq
i}K_{ik}\gamma_{ik}\geq\displaystyle\sum_{k\neq
i}K_{ik}\gamma_{ik}=\displaystyle\sum_{k\neq
i}\left|K_{ik}\gamma_{ik}\right|=\displaystyle\sum_{j\neq
i}\left|\mathcal{J}_{ij}\right|.$ (6)

According to the Gershgorin circle theorem, every eigenvalue of the Jacobian
matrix lies within at least one of the $N$ Gershgorin discs
$D_{i}(\mathcal{J}_{ii},r_{i})=\left\\{z\in\mathbb{C}|\,\,\left|z-\mathcal{J}_{ii}\right|\leq
r_{i}\right\\}$ with $r_{i}=\sum_{k\neq i}|\mathcal{J}_{ki}|$ in the complex
plane. Relation (6) indicates that all $N$ Gershgorin discs lie in the left
half of the complex plane, i.e. in
$\left\\{z\in\mathbb{C}\left|\,\mbox{Re}\left(z\right)\leq 0\right.\right\\}$,
because the centers of the discs $(\mathcal{J}_{ii},0)$ lie on the negative
real axis and the radius of the discs $r_{i}=\sum_{k\neq
i}|\mathcal{J}_{ki}|\leq|\mathcal{J}_{ii}|$. The discs touch the imaginary
axis from the half plane
$\left\\{z\in\mathbb{C}\left|\,\mbox{Re}\left(z\right)\leq 0\right.\right\\}$
only if $\beta_{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 $\bm{x}=\bm{x}^{*}$.

######  Remark 2.4 (Homogeneous nodal dynamics and graph Laplacian)

We remark that in case of identical intrinsic nodal dynamics at all nodes, a
weighted Laplacian matrix $\mathcal{L}$ of the interaction graph explicitly
enters the linearized response dynamics of the network (3). Assuming
$\beta_{i}=\beta\in\mathbb{R}$ for all nodes $i$, we can express the Jacobian
matrix as

$\mathcal{J}=-\beta\bm{1}-\mathcal{L},$ (7)

where the weighted graph Laplacian $\mathcal{L}$ is defined as

$\mathcal{L}_{ij}:=\left\\{\begin{array}[]{ll}\displaystyle\sum_{k\neq
i}K_{ik}\gamma_{ik}&\mathrm{for}\quad i=j\\\
-K_{ij}\gamma_{ij}&\mathrm{for}\quad i\neq j.\end{array}\right.$ (8)

Here $K_{ij}\gamma_{ij}$ is considered a weight of edge $(i,j)$, containing
the coupling strength $K_{ij}$ and the linearized coupling function
$\gamma_{ij}$ at the fixed point.

######  Remark 2.5 (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 $\gamma_{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
$\gamma_{ji}$ for the edge $(j,i)$ with the opposite direction. This asymmetry
leads to an asymmetric weighted graph Laplacian (8). However, for undirected
networks with symmetric strengths ($K_{ij}=K_{ji}$) and symmetric
sensitivities of coupling functions ($\gamma_{ij}=\gamma_{ji}$), or more
generally, a symmetric combination $K_{ij}\gamma_{ij}=K_{ji}\gamma_{ji}$, the
weighted graph Laplacian $\mathcal{L}$ is symmetric. Its eigenvectors thus
form an orthogonal basis which allows us to solve for the linear network
responses in (3) and (7) by expressing the response vector $\bm{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 $\varepsilon>0$ and frequency $\omega>0$ at node $k$, i.e.
$D^{(k)}_{i}(t)=\delta_{ik}\varepsilon e^{\imath(\omega t+\varphi)}$, where
$\delta_{ik}$ is the Kronecker delta with $\delta_{ik}=1$ if $i=k$ and
$\delta_{ik}=0$ otherwise. We solve for the linear network response vector
$\bm{X}^{(k)}(\omega,t)$ governed by

$\dot{\bm{X}}^{(k)}=-\beta\bm{X}^{(k)}-\mathcal{L}\bm{X}^{(k)}+\bm{D}^{(k)}(t).$
(9)

Expressing the response vector in the constant eigenbasis of Laplacian

$\bm{X}^{(k)}(t)=\sum_{\ell=0}^{N-1}c^{[\ell]}(t)\bm{v}^{[\ell]}$ (10)

and the time-dependent coefficients $c^{[\ell]}(t)$ and exploiting
orthogonality of the Laplacian eigenvectors, we obtain the ordinary
differential equations

$\dot{c}^{[\ell]}=(-\beta-\lambda^{[\ell]})c^{[\ell]}+\varepsilon
v^{[\ell]}_{k}e^{\imath(\omega t+\varphi)}$ (11)

for the coefficients. Here the $N$ Laplacian eigenvalues $\lambda^{[\ell]}$
and eigenvectors $\bm{v}^{[\ell]}$ are indexed according to
$0=\lambda^{[0]}\leq\lambda^{[1]}\leq\cdots\leq\lambda^{[N-1]}$. 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 (11) for the coefficients is

$c^{[\ell]}(\omega,t)=\left(\bm{X}_{0}^{(k)}\cdot\bm{v}^{[\ell]}\right)e^{\left(-\beta-\lambda^{[\ell]}\right)t}+\dfrac{\varepsilon
v^{[\ell]}_{k}e^{\imath\varphi}}{\beta+\lambda^{[\ell]}+\imath\omega}\left(-e^{\left(-\beta-\lambda^{[\ell]}\right)t}+e^{\imath\omega
t}\right),$ (12)

where $\bm{X}_{0}^{(k)}$ denotes the initial response vector at $t=0$ given
node $k$ is perturbed. The linear network response is thus given by (10) and
(12).

######  Remark 2.6 (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
$\bm{D}(t)$ are arbitrary time series, the linear network response can be
obtained by summing up the single-signal single-frequency response (10 and 12)
over all frequency components $\omega$ and all perturbation signals $k$ by
resorting to the linearity of the dynamics (3)

$\bm{X}(t)=\sum_{k}\left(\bm{X}^{(k)}_{\omega=0}(t)+\sum_{\omega>0}\bm{X}^{(k)}(\omega,t)\,\mathrm{d}\omega\right).$
(13)

Here $\bm{X}^{(k)}_{\omega=0}(t)$ denotes the linear response to a constant
($\omega=0$) perturbation $\varepsilon$, which shares the same form as
$\bm{X}^{(k)}(\omega,t)$ with $\omega=0$ for $\beta>0$ but not for $\beta=0$.
Moreover, the sum becomes an integral for continuously distributed
frequencies, with $\bm{X}^{(k)}(\omega,t)$ replaced by the associated response
density per unit frequency.

For systems where nodal damping vanishes ($\beta=0$), the coefficient for the
$0$-th eigenmode $c^{[0]}$ diverges for a constant perturbation as
$\omega\rightarrow 0$, so it needs to be solved separately and takes the form
of $\bm{X}_{0}^{(k)}\cdot\bm{v}^{[0]}+\varepsilon v_{k}^{[0]}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.7 (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 (7)
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.5. If the higher-order time derivatives of the state
variables has homogeneous coefficients for individual nodes, i.e. the response
dynamics to a perturbation $\bm{D}(t)$ has the form of

$\sum_{d}\kappa_{d}D_{t}^{d}\bm{X}=-\mathcal{L}\bm{X}+\bm{D}(t),$ (14)

an explicit solution of the linear responses in the eigenbasis of
$\mathcal{L}$ can still be obtained following the routine in Remark 2.5, if
the corresponding ODEs for the time-dependent coefficients

$\sum_{d}\kappa_{d}D_{t}^{d}{c}^{[\ell]}=-\lambda^{[\ell]}c^{[\ell]}+\bm{D}(t)\cdot\bm{v}^{[\ell]}$
(15)

are solvable. Here $\kappa_{d}\in\mathbb{R}$ are constant coefficients and we
use Euler’s notation for derivatives, where $D_{t}^{d}x$ denotes the $d$-th
time derivative of variable $x$.

In summary, if a networked dynamical system consisting of $N$ units

  1. 1.

has a fixed point at $\bm{x}=\bm{x}^{*}$,

  2. 2.

in the neighbourhood of $\bm{x}=\bm{x}^{*}$ the intrinsic nodal dynamics
$h_{i}(x_{i})$ gives homogeneous non-positive feedback to the respective
nodes, i.e.
$\beta_{i}:=-\left.\frac{\mathrm{d}h_{i}}{\mathrm{d}x_{i}}\right|_{\bm{x}=\bm{x}^{*}}=\beta\geq
0$ for all $i$, and

  3. 3.

the coupling function $g_{ij}(x_{j}-x_{i})$ is diffusive and the coupling
term’s sensitivity to small changes at $\bm{x}=\bm{x}^{*}$ is symmetric and
non-negative, i.e. $K_{ij}=K_{ji}$ and
$\gamma_{ij}:=\left.\frac{\mathrm{d}g_{ij}}{\mathrm{d}(x_{j}-x_{i})}\right|_{\bm{x}=\bm{x}^{*}}=\gamma_{ji}\geq
0$ for all edges $(i,j)$,

then i) the dynamical system is at least neutrally stable (2.2) and ii) a
symmetric weighted graph Laplacian (8) arises in the network response dynamics
(9) and enables the expression of linear network responses in the Laplacian
eigenbasis (10) and (12) (Remark 2.4 and Remark 2.5). 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.

##  3 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.

###  3.1 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 model111The DC power flow model here
must not be confused with models for high-voltage direct-current (HVDC)
transmission grids, which actually uses DC, as opposed to alternating current
(AC), for the transmission of electrical power. 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

$S_{ij}=U_{j}I^{*}_{ij}=U_{j}\left(\dfrac{U_{j}-U_{i}}{Z_{ij}}\right)^{*}=\dfrac{|U_{j}|e^{\imath\theta_{j}}(|U_{j}|e^{-\imath\theta_{j}}-|U_{i}|e^{-\imath\theta_{i}})}{R_{ij}-\imath
X_{ij}},$ (16)

where $U_{j}=|U_{j}|e^{\imath\theta_{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}+\imath
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}+\imath 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: [Perfect voltage
support] The voltage amplitude is constant and identical for each node in the
power grid network, $|U_{i}|\equiv U$ for all $i$, and the management of the
reactive power is not considered. [Lossless lines] Transmission losses on the
lines are neglected, implying that line resistances are negligible compared to
line reactances: $R_{ij}/X_{ij}\rightarrow 0$ for all lines $(i,j)$. [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
$\pi/2$ such that $\sin(\theta_{j}-\theta_{i})\approx\theta_{j}-\theta_{i}$
and $\cos(\theta_{j}-\theta_{i})\approx 1$.

With the above mentioned assumptions in mind, the complex power flow (16)
between neighboring nodes simplifies to

$S_{ij}=\dfrac{U^{2}(1-e^{\imath(\theta_{j}-\theta_{i})})}{-\imath
X_{ij}}=\dfrac{U^{2}}{X_{ij}}\left(\theta_{j}-\theta_{i}\right)=P_{ij}.$ (17)

Here the complex power flow $S_{ij}$ naturally reduces to the active power
flow $P_{ij}$ since the imaginary part vanishes. The equation (17) resembles
the expression of the direct current carried by a “resistor” $X_{ij}/U^{2}$
influenced by a “voltage drop” $\theta_{j}-\theta_{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

$P_{i}=\sum_{j=1}^{N}P_{ij}=\sum_{j=1}^{N}K_{ij}\left(\theta_{i}-\theta_{j}\right)$
(18)

over all connected units. Equation (18) is the core of the DC power flow model
as it yields the pattern of power flows $K_{ij}\theta_{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,
$\bm{P}:=(P_{1},\cdots,P_{N})$ and
$\bm{\theta}:=(\theta_{1},\cdots,\theta_{N})$, we express the linear relation
between them by a weighted graph Laplacian $\mathcal{L}$ introduced in Section
2.2, i.e.

$\bm{P}=\mathcal{L}\bm{\theta},$ (19)

Here $\mathcal{L}$ is defined similarly as in (8), only with
$\gamma_{ij}\equiv 1$ for all edges $(i,j)$.

Assuming that the power transmission network runs at a normal operation state
where the voltage angles $\bm{\theta}^{*}$ are stationary at all nodes, the
fixed voltage angle differences determine a specific power flow pattern
$\bm{P}^{*}$ across the network through the linear operator $\mathcal{L}$ such
that $\bm{P}^{*}=\mathcal{L}\bm{\theta}^{*}$. If the nodal power injections
are perturbed as $\bm{P}(t)=\bm{P}^{*}+\bm{D}(t)$ by a shift vector
$\bm{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 $\bm{\theta}=\bm{\theta}^{*}+\bm{\Theta}$ change accordingly
through $\mathcal{L}$ due to the linear operation in (19). The response vector
of the voltage angles $\bm{\Theta}$ is given by

$\bm{\Theta}(t)=\mathcal{L}^{+}\bm{D}(t),$ (20)

where $\mathcal{L}^{+}$ denotes the Moore-Penrose inverse of the weighted
graph Laplacian $\mathcal{L}$.

###### Remark 3.1

The weighted graph Laplacian $\mathcal{L}$ defined in (8) is singular since it
always has an eigenvalue $\lambda_{0}=0$ with the corresponding eigenvector
$\bm{v}_{0}=(1,\cdots,1)$ satisfying $\mathcal{L}\bm{v}_{0}=\bm{0}$ by
construction. Therefore, to compute the voltage angle shifts (20) we need the
generalized inverse matrix $\mathcal{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 $\theta_{k}$ of one of the units $k$ as
the reference for voltage angle, i.e. by setting
$\theta_{i}\rightarrow\theta_{i}-\theta_{k}$ for all $i$. If the network is
connected, then the $(N-1)$-dimensional Laplacian matrix is invertible and
both of the matrix equations, (19) and (20), have a unique solution
respectively.

###  3.2 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 $\Omega_{0}^{m}$ corresponding to the
nominal grid frequency $\Omega_{0}=2\pi\times 50$ or $2\pi\times 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)]:

$I\ddot{\theta}^{\text{m}}+D^{\text{m}}\dot{\theta}^{\text{m}}=T^{\text{m}}-T^{\text{e}},$
(21)

where $\theta^{\text{m}}$ denotes the mechanical rotor angle deviation from
the rotating reference frame $\Omega_{0}t$, $I$ denotes the moment of inertia
of the rotor and the connected turbine, $D^{\text{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^{\text{m}}$ and $T^{\text{e}}$ denote respectively the net mechanical
torque and the counteracting electromagnetic torque acting on the rotor.

For deviations of the angular velocity $\dot{\theta}^{\text{m}}$, the local
frequency deviation is small compared to the nominal grid frequency
$\Omega_{0}^{\text{m}}$, i.e.
$(\dot{\theta}^{\text{m}}+\Omega_{0}^{\text{m}})^{-1}\approx(\Omega_{0}^{\text{m}})^{-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=(\dot{\theta}^{\text{m}}+\Omega_{0}^{\text{m}})^{-1}P\approx(\Omega_{0}^{\text{m}})^{-1}P$.
and using the relations Also introducing the effective quantities
$\theta(t)=\theta^{\text{m}}(t)/(p/2)$ and
$\Omega_{0}=\Omega^{\text{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:

$M\ddot{\theta}+\tilde{D}\dot{\theta}=P^{\text{m}}-P^{\text{e}}.$ (22)

Here $M:=I\Omega_{0}/(p/2)^{2}$ and
$\tilde{D}:=D^{\text{m}}\Omega_{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 (22), $P^{\text{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^{\text{e}}$ denotes the electrical power injected to the grid by the
synchronous machine. In the oscillator model we again assume perfect voltage
support (subsection 3.1) and lossless transmission lines (subsection 3.1),
which yield $P^{\text{e}}=\sum_{j=1}^{N}K_{ij}(\theta_{j}-\theta_{i})$ [cf.
$P_{i}$ in the DC power flow model (18)]. Putting everything together, we
obtain the governing equations of the oscillator model of AC power grids

$\ddot{\theta}_{i}=P_{i}-\alpha_{i}\dot{\theta}_{i}+\sum_{j=1}^{N}K_{ij}\sin(\theta_{j}-\theta_{i}),\quad\text{for
}i\in\\{1,\cdots,N\\}$ (23)

with the parameters $P_{i}:=P_{i}^{\text{m}}/M_{i}$,
$\alpha_{i}:=\tilde{D}_{i}/M_{i}$, and $K_{ij}:=U^{2}/X_{ij}$.

######  Proposition 3.2 (Linear stability of the oscillator model)

AC power grid systems described by the oscillator model (23) with underlying
interaction topology $G(V,E)$ is at least neutrally stable at a fixed point
$\bm{\theta}^{*}$, if all edges are not overloaded, i.e.
$|\theta_{j}^{*}-\theta_{i}^{*}|\leq\frac{\pi}{2}$ for all $(i,j)\in E$.

###### Proof 3.3

At a fixed point of the system $\bm{\theta}=\bm{\theta}^{*}$, a small
deviation of the oscillators’ angles
$\bm{\Theta}:=\bm{\theta}-\bm{\theta}^{*}$ follows the linear dynamics

$\displaystyle\dfrac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix}\bm{\Theta}\\\
\dot{\bm{\Theta}}\end{pmatrix}=\mathcal{J}\begin{pmatrix}\bm{\Theta}\\\
\dot{\bm{\Theta}}\end{pmatrix},$ (24)

where the Jacobian matrix $\mathcal{J}\in\mathbb{R}^{2N\times 2N}$ of the
$2N$-dimensional dynamical system is given by

$\displaystyle\mathcal{J}=\begin{pmatrix}\bm{0}_{N}&\bm{I}_{N}\\\
-\mathcal{L}&-\mathcal{A}.\end{pmatrix}$ (25)

Here $\mathcal{L}$ is a weighted graph Laplacian as defined in (8) with
$\gamma_{ij}=\cos(\theta^{*}_{j}-\theta^{*}_{i})$, $\mathcal{A}$ is an
$N\times N$ diagonal matrix with $\mathcal{A}_{ii}:=\alpha_{i}$, and
$\bm{0}_{N}$ and $\bm{I}_{N}$ are respectively the $N\times N$ zero matrix and
identity matrix.

Let $\bm{w}=(\bm{w}_{1},\bm{w}_{2})\in\mathbb{C}^{2N}$ with
$\bm{w}_{1},\bm{w}_{2}\in\mathbb{C}^{N}$ be an eigenvector of $\mathcal{J}$
corresponding to eigenvalue $\mu\in\mathbb{C}$. By definition we have
$\mathcal{J}\bm{w}=\mu\bm{w}$, which by writing (24) as a second order
differential equation implies

$\mu^{2}\bm{w}_{1}+\mu\mathcal{A}\bm{w}_{1}+\mathcal{L}\bm{w}_{1}=\bm{0}.$
(26)

Multiplying both sides of the equation above with the conjugate transpose,
$\bm{w}_{1}^{\dagger}$ of $\bm{w}_{1}$ from the left, we obtain an expression
of the eigenvalue
$\mu=\left(-\chi_{2}\pm\sqrt{\chi_{2}^{2}-4\chi_{1}\chi_{3}}\right)/2\chi_{1}$,
with $\chi_{1}=\bm{w}^{\dagger}_{1}\bm{w_{1}}\geq 0$,
$\chi_{2}=\bm{w}^{\dagger}_{1}\mathcal{A}\bm{w}_{1}\geq 0$ and
$\chi_{3}=\bm{w}^{\dagger}_{1}\mathcal{L}\bm{w}_{1}\geq 0$. $\chi_{2}$ and
$\chi_{3}$ are non-negative since $\mathcal{A}_{ii}=\alpha_{i}>0$ and
$\mathcal{L}$ is positive semi-definite because $\gamma_{ij}\geq 0$ is ensured
through $|\theta_{j}^{*}-\theta_{i}^{*}|\leq\frac{\pi}{2}$ for all $(i,j)\in
E$ (see 2.2 and Remark 2.4) [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.4 (Neutral stability and global phase shift)

Connected AC power grids described by the oscillator model (23) is neutrally
stable at a fixed point only when the deviation $\bm{\Theta}$ is a global
phase shift. Because (26) indicates that the Jacobian eigenvalue $\mu=0$ only
when $\chi_{3}=0$, which implies $\mathcal{L}\bm{w}_{1}=\bm{0}$ and thus
$\bm{w}_{1}$ lies in the Laplacian eigenspace corresponding to the eigenvalue
$\lambda^{[0]}=0$. For connected graphs there is only one eigenvector
$\bm{v}^{[0]}=(1,\cdots,1)$ corresponding to $\lambda^{[0]}=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.5 (LRT of the oscillator model and homogeneous nodal
damping)

Consider an AC power grid oscillator model with arbitrary topology (23) with
homogeneous nodal damping $\alpha_{i}=\alpha\geq 0$ for all nodes $i$. Then
the network-wide linear responses to arbitrary external perturbations near a
normal operation state $\bm{\theta}^{*}$ can be expressed explicitly in the
eigenbasis of a weighted graph Laplacian: i) The network response to time-
independent distributed perturbations $\bm{D}^{*}$ is

$\bm{\Theta}(t)=\bm{D}^{*}\cdot\bm{v}^{[0]}\left(\dfrac{e^{-\alpha
t}}{\alpha^{2}}-\dfrac{1}{\alpha^{2}}+\dfrac{t}{\alpha}\right)\bm{v}^{[0]}+\sum_{\ell=1}^{N-1}\dfrac{\bm{D}^{*}\cdot\bm{v}^{[\ell]}}{\lambda^{[\ell]}}\left(\dfrac{\Delta_{-}^{[\ell]}e^{\Delta_{+}^{[\ell]}t}-\Delta_{+}^{[\ell]}e^{\Delta_{-}^{[\ell]}t}}{2\eta^{[\ell]}}+1\right)\bm{v}^{[\ell]},$
(27)

and ii) the network response to a single sinusoidal perturbation given by
$\bm{D}(t)$ with $D_{i}(t)=\delta_{ik}\varepsilon e^{\imath(\omega
t+\varphi)}$ is

$\bm{\Theta}^{(k)}(t)=\sum_{\ell=0}^{N-1}\dfrac{\varepsilon
v_{k}^{[\ell]}e^{\imath\varphi}}{-\omega^{2}+\imath\alpha\omega+\lambda^{[\ell]}}\left[\dfrac{\left(\Delta^{[\ell]}_{-}-\imath\omega\right)e^{\Delta^{[\ell]}_{+}t}-\left(\Delta^{[\ell]}_{+}-\imath\omega\right)e^{\Delta^{[\ell]}_{-}t}}{2\eta^{[\ell]}}+e^{\imath\omega
t}\right]\bm{v}^{[\ell]}$ (28)

with $\Delta^{[\ell]}_{\pm}:=-\alpha/2\pm\eta^{[\ell]}$ and
$\eta^{[\ell]}:=\sqrt{\alpha^{2}/4-\lambda^{[\ell]}}$.

###### Proof 3.6

Vectorizing the linear response $\bm{\Theta}(t)$ of the system (23) to a
perturbation vector $\bm{D}(t)$, we obtain the matrix equation describing the
response dynamics of the oscillator model of AC power grids

$\ddot{\bm{\Theta}}+\bm{\alpha}\circ\dot{\bm{\Theta}}=-\mathcal{L}\bm{\Theta}+\bm{D}(t),$
(29)

where $\bm{\alpha}:=(\alpha_{1},\cdots,\alpha_{N})$ denotes the vector of
damping parameters and $``\circ"$ denotes the Schur (element-wise) product of
two vectors. Let $\alpha_{i}=\alpha$ for all $i$, the term
$\bm{\alpha}\circ\dot{\bm{\Theta}}$ reduces to a scalar multiplication
$\alpha\dot{\bm{\Theta}}$, thereby all terms involving the variable
$\bm{\Theta}$ in equation (29) can be expressed as linear combinations of
Laplacian eigenvectors. Because the $\mathcal{L}$ here is real and symmetric
so that we can write
$\bm{\Theta}(t)=\sum_{\ell=0}^{N-1}c^{[\ell]}(t)\bm{v}^{[\ell]}$. Using the
same trick as in Remark 2.5, we obtain equations for the time-dependent
coefficients $c^{[\ell]}(t)$ given by

$\ddot{c}^{[\ell]}+\alpha\dot{c}^{[\ell]}+\lambda^{[\ell]}c^{[\ell]}=\bm{D}(t)\cdot\bm{v}^{[\ell]}\quad\text{for
}\ell\in\\{0,\cdots,N-1\\}.$ (30)

Assuming at $t=0$ the AC power grid system operates normally at the fixed
point $\bm{\theta}^{*}$, we have initial conditions for the coefficients
$c^{[\ell]}(0)=0$ and $\dot{c}^{[\ell]}(0)=0$ for all $\ell$, thereby obtain
explicitly solutions for the coefficients and for the network linear
responses.

For perturbations independent of time, $\bm{D}(t)=\bm{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 (30):

$\bm{\Theta}(t)=\bm{D}^{*}\cdot\bm{v}^{[0]}\left(\dfrac{e^{-\alpha
t}}{\alpha^{2}}-\dfrac{1}{\alpha^{2}}+\dfrac{t}{\alpha}\right)\bm{v}^{[0]}+\sum_{\ell=1}^{N-1}\dfrac{\bm{D}^{*}\cdot\bm{v}^{[\ell]}}{\lambda^{[\ell]}}\left(\dfrac{\Delta_{-}^{[\ell]}e^{\Delta_{+}^{[\ell]}t}-\Delta_{+}^{[\ell]}e^{\Delta_{-}^{[\ell]}t}}{2\eta^{[\ell]}}+1\right)\bm{v}^{[\ell]},$
(31)

with $\Delta^{[\ell]}_{\pm}:=-\alpha/2\pm\eta^{[\ell]}$ and
$\eta^{[\ell]}:=\sqrt{\alpha^{2}/4-\lambda^{[\ell]}}$. For time-dependent
perturbations $\bm{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.6. Similarly, we let $D_{i}(t)=\delta_{ik}\varepsilon
e^{\imath(\omega t+\varphi)}$ and obtain the oscillator model’s linear
response to a sinusoidal signal at node $k$ as

$\bm{\Theta}^{(k)}(t)=\sum_{\ell=0}^{N-1}\dfrac{\varepsilon
v_{k}^{[\ell]}e^{\imath\varphi}}{-\omega^{2}+\imath\alpha\omega+\lambda^{[\ell]}}\left[\dfrac{\left(\Delta^{[\ell]}_{-}-\imath\omega\right)e^{\Delta^{[\ell]}_{+}t}-\left(\Delta^{[\ell]}_{+}-\imath\omega\right)e^{\Delta^{[\ell]}_{-}t}}{2\eta^{[\ell]}}+e^{\imath\omega
t}\right]\bm{v}^{[\ell]}.$ (32)

######  Remark 3.7 (Low dissipation regime and grid eigenfrequencies)

In case the dissipation in the system is low enough such that
$\alpha<2\sqrt{\lambda^{[\ell]}}$ for the $\ell$-th eigenvalue, in the
solution of the linear response (28) the corresponding $\eta^{[\ell]}$ 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^{-\frac{\alpha}{2}t}$, i.e., with a time constant
$\tau=2/\alpha$. Such intrinsic oscillation modes can also be identified by
looking at the eigenvalues of the Jacobian matrix at the fixed point. Let
$\bm{w}_{1}$ in (26) be the $\ell$-th eigenvector of the Laplacian
$\mathcal{L}$, we can see the corresponding Jacobian eigenvalue
$\mu^{[\ell]}=-\alpha/2\pm\sqrt{\alpha^{2}/4-\lambda^{[\ell]}}=\Delta^{[\ell]}_{\pm}$,
which indicates the corresponding eigenfrequency
$\omega^{[\ell]}_{\text{eigen}}:=\mbox{Im}[\mu^{[\ell]}]=\sqrt{\lambda^{[\ell]}-\alpha^{2}/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
$\alpha<2\sqrt{\lambda^{[\ell]}}$ for all $N-1$ Laplacian eigenvalues.

##  4 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.

###  4.1 Frequency regimes of steady-state response patterns

After a transient phase characterized by a dissipation-related time constant
$\tau=2/\alpha$ (cf. Remark 3.7), the perturbed power grid systems reside in a
second regime of network responses, where the entire network respond
periodically to perturbation signals for $t\gg\tau$ (see 27 and 28). 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 (23), the steady-state responses
to time-independent perturbations $\bm{D}(t)=\bm{D}^{*}$, i.e. with a
perturbation frequency $\omega=0$, near a normal operation state
$\bm{\theta}^{*}$ are constituted by a homogeneous shift of grid frequency

$\delta\dot{\theta}_{i}=\dfrac{1}{N\alpha}\displaystyle\sum_{j=1}^{N}D^{*}_{j}\quad\text{for
}i\in\\{1,\cdots,N\\},$ (33)

and a topology-dependent phase shift

$\delta\bm{\theta}=-\dfrac{1}{N\alpha^{2}}\displaystyle\sum_{i=1}^{N}D^{*}_{i}\bm{1}+\sum_{\ell=1}^{N-1}\dfrac{\bm{D}^{*}\cdot\bm{v}^{[\ell]}}{\lambda^{[\ell]}}\bm{v}^{[\ell]}.$
(34)

###### Proof 4.2

By definition, the steady-state response patterns become clear by
investigating the asymptotic behaviour of the responses as
$t\rightarrow\infty$. For responses to a constant perturbation vector
$\bm{D}(t)=\bm{D}^{*}$ (27) the steady-state response reads

$\bm{\Theta}(t)\overset{t\rightarrow\infty}{\sim}\bm{D}^{*}\cdot\bm{v}^{[0]}\left(-\dfrac{1}{\alpha^{2}}+\dfrac{t}{\alpha}\right)\bm{v}^{[0]}+\sum_{\ell=1}^{N-1}\dfrac{\bm{D}^{*}\cdot\bm{v}^{[\ell]}}{\lambda^{[\ell]}}\bm{v}^{[\ell]},$
(35)

which consists of two characteristic patterns: i) the phases of all units
drift away from the normal operation state with a constant angular velocity
$(\bm{D}^{*}\cdot\bm{v}^{[0]})v^{[0]}_{i}/\alpha$, and ii) a time-independent
and unit-specific phase shift. Since $\bm{v}^{[0]}=\frac{1}{\sqrt{N}}\bm{1}$,
the former pattern represents the constant homogeneous shift of grid frequency
(33) and the latter the topology-dependent phase shift (34).

###### Remark 4.3

The global grid frequency shift (33) is a consequence of the imbalance between
power injected into and drawn from the power grid system, which is imposed by
the constant perturbation $\bm{D}^{*}$. The topology-dependent phase shifts
(34) 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(\theta_{j}-\theta_{i})$ on line $(i,j)$, $\delta
L_{ij}=\cos(\theta_{j}^{*}-\theta_{i}^{*})(\delta\theta_{j}-\delta\theta_{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.4 (Steady-state response pattern for a sinusoidal
perturbation)

For AC power grids with arbitrary topologies (23), the steady-state responses
to a sinusoidal perturbation at node $k$, i.e. $\bm{D}(t)$ with
$D_{i}(t)=\delta_{ik}\varepsilon e^{\imath(\omega t+\varphi)}$, near a normal
operation state $\bm{\theta}^{*}$ are constituted by a homogeneous phase shift

$\delta\theta_{i}=\dfrac{\imath\varepsilon e^{\imath\varphi}}{\alpha\omega
N}\quad\text{for }i\in\\{1,\cdots,N\\},$ (36)

and sinusoidal responses at each node with the same frequency $\omega$ and a
characteristic complex amplitude

$R_{i}^{(k)}(\omega):=\sum_{\ell=0}^{N-1}\dfrac{v_{k}^{[\ell]}v_{i}^{[\ell]}}{-\omega^{2}+\imath\alpha\omega+\lambda^{[\ell]}}.$
(37)

for node $i$.

###### Proof 4.5

The steady-state response to a sinusoidal perturbation at a single node $k$ is
obtained by studying the asymptotic behaviour of the responses (28) has the
form of

$\bm{\Theta}^{(k)}(t)\overset{t\rightarrow\infty}{\sim}\dfrac{\imath\varepsilon
e^{\imath\varphi}}{\alpha\omega N}\bm{1}+e^{\imath(\omega
t+\varphi)}\sum_{\ell=0}^{N-1}\dfrac{\varepsilon
v_{k}^{[\ell]}}{-\omega^{2}+\imath\alpha\omega+\lambda^{[\ell]}}\bm{v}^{[\ell]},$
(38)

which is composed of a homogeneous phase shift $\frac{\imath\varepsilon
e^{\imath\varphi}}{\alpha\omega N}$ and a driven oscillation at each node.
Each node’s angular variable $\theta_{i}$ changes at the same frequency as the
perturbation frequency $\omega$, but with a complex amplitude

$R_{i}^{(k)}(\omega)=\sum_{\ell=0}^{N-1}\dfrac{v_{k}^{[\ell]}v_{i}^{[\ell]}}{-\omega^{2}+\imath\alpha\omega+\lambda^{[\ell]}}.$
(39)

######  Remark 4.6 (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_{i}^{(k)}$ 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_{i}^{(k)}$ determines the maximal deviation of the local
grid frequency at node $i$ caused by a perturbation at node $k$ through
$|\delta\dot{\theta}_{i}|=|\dot{\Theta}^{(k)}_{i}|=\varepsilon\omega|R_{i}^{(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

$A_{i}^{(k)}:=\omega\left|R_{i}^{(k)}\right|$ (40)

for each node $i$ and discuss in detail its distinctive spatial distributions
in different frequency regimes.

Figure 1: Steady-state response patterns exhibit three frequency regimes. (a)
Relative response strength $A_{i}^{*(k)}$ (see Remark 4.10 for definition) of
all $80$ nodes in an example power grid network across all three frequency
regimes (homogeneous bulk, resonant, and localized responses). Vertical grey
lines represent the $N-1$ eigenfrequencies. (b1,c1,d1) Qualitatively different
dependencies of $A_{i}^{*(k)}$ on the graph-theoretic distance $d:=d(k,i)$
with three representative driving frequencies $\omega/2\pi\in\\{0.1,2,10\\}$
Hz of three frequency regimes. The exponential dependence of $A_{i}^{*(k)}$ on
$d$ is illuatrated in the inset of d1. (b2,c2,d2) Distinctive response
patterns for the three driving frequencies, corresponding to (b1,c1,d1). The
curves in (a) are color-coded with the distance $d$, and the discs in
(b1-b2,c1-c2,d1-d2) are color-coded with the relative response strength
$A_{i}^{*(k)}$. The black square marks the perturbed node. Network settings
are the same as Fig. 2 in Zhang et al. (2019).

####  4.1.1 Regime of grid resonances

As suggested in Remark 3.7, 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
$\lambda^{[\ell]}>\alpha^{2}/4$, the oscillatory power grid system resonates
if the perturbation frequency matches the corresponding eigenfrequency
$\omega_{\text{eigen}}^{[\ell]}=\sqrt{\lambda^{[\ell]}-\alpha^{2}/4}$. Driven
at frequencies close to an eigenfrequency $\omega_{\text{eigen}}^{[\ell]}$,
the network responses exhibit large amplitudes since the rationalized
denominator $(\omega^{2}-\lambda^{[\ell]})^{2}-\alpha^{2}\omega^{2}$ for the
corresponding $\ell$-th oscillation mode is minimized. However, the response
amplitudes vary greatly for different nodes in the network due to the factor
$v_{k}^{[\ell]}v_{i}^{[\ell]}$. 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.7).

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).

####  4.1.2 Homogeneous responses: the low-frequency regime

The network response pattern for lower perturbation frequencies, i.e. the ones
lower than the smallest eigenfrequency
$\omega_{\text{eigen}}^{[1]}=\sqrt{\lambda^{[1]}-\alpha^{2}/4}$, can be
understood by investigating the asymptotic behaviour of the nodal response
strength $A_{i}^{(k)}=\omega|R_{i}^{(k)}|$ as the $\omega\rightarrow 0$.

######  Proposition 4.7 (Homogeneous responses at the low frequency limit)

As the perturbation frequency $\omega\rightarrow 0$, the steady-state response
strength at each node of an AC power grid system with an arbitrary topology
(23) approaches a constant value, i.e.

$A_{i}^{(k)}\overset{\omega\rightarrow 0}{\sim}\frac{1}{N\alpha}.$ (41)

###### Proof 4.8

Considering (37) and $\lambda^{[0]}=0$, we derive the asymptotic behaviour of
the real part and the imaginary part of the nodal response factor
$R_{i}^{(k)}$ as $\omega\rightarrow 0$:

$\displaystyle\mbox{Re}\left[R_{i}^{(k)}\right]$
$\displaystyle=\sum_{\ell=0}^{N-1}\frac{v^{[\ell]}_{k}v^{[\ell]}_{i}\left(-\omega^{2}+\lambda^{[\ell]}\right)}{\left(-\omega^{2}+\lambda^{[\ell]}\right)^{2}+\alpha^{2}\omega^{2}}\overset{\omega\rightarrow
0}{\sim}-\dfrac{1}{N\alpha^{2}}+\displaystyle\sum_{\ell=1}^{N-1}\dfrac{v^{[\ell]}_{k}v^{[\ell]}_{i}}{\lambda^{[\ell]}}$
(42) $\displaystyle\mbox{Im}\left[R_{i}^{(k)}\right]$
$\displaystyle=\sum_{\ell=0}^{N-1}\frac{v^{[\ell]}_{k}v^{[\ell]}_{i}\left(-\alpha\omega\right)}{\left(-\omega^{2}+\lambda^{[\ell]}\right)^{2}+\alpha^{2}\omega^{2}}\overset{\omega\rightarrow
0}{\sim}-\dfrac{1}{N\alpha\omega}.$ (43)

As the asymptotic behaviour of the response strength
$A_{i}^{(k)}=\omega|R_{i}^{(k)}|$ is dominated by the imaginary part, we have

$A_{i}^{(k)}=\omega\left|R_{i}^{(k)}\right|\overset{\omega\rightarrow
0}{\sim}\omega\cdot\dfrac{1}{N\alpha\omega}=\dfrac{1}{N\alpha}.$ (44)

######  Remark 4.9 (Consistency with the homogeneous grid frequency shift at
$\omega=0$)

The homogeneous response strength at each node as $\omega\rightarrow 0$
suggests a global shift of grid frequency inversely proportional to the
network size $N$ and the system dissipation parameter $\alpha$. This result is
quantitatively consistent with the homogeneous grid frequency shift induced by
constant perturbations, as discussed in 4.1 with $D_{i}^{*}=\delta_{ik}$ if
node $k$ is the perturbed one.

######  Remark 4.10 (Relative nodal response strength)

The homogeneous nodal response at the low frequency limit (41) serves as a
reference value for the nodal response strengths of a sinusoidally driven
network. Therefore, by normalizing $A_{i}^{(k)}$ with its low frequency limit
$\frac{1}{N\alpha}$, we define the relative nodal response strength

$\displaystyle
A_{i}^{*(k)}:=\dfrac{A_{i}^{(k)}}{\displaystyle\lim_{\omega\rightarrow
0}A_{i}^{(k)}}=N\alpha A_{i}^{(k)}$ (45)

so that the nodal response strengths can be compared across networks with
different sizes and dissipation values.

####  4.1.3 Localized responses: the high-frequency regime

Similarly, we investigate the network response pattern in the high-frequency
regime where $\omega>\omega_{\text{eigen}^{[N-1]}}$ by observing the
asymptotic behaviour of the nodal response strengths $A_{i}^{(k)}$ as the
perturbation frequency becomes sufficiently large and approaches infinity.

######  Proposition 4.11 (Localized response patterns in the high-frequency
regime)

As the perturbation frequency $\omega\rightarrow\infty$, the steady-state
nodal response strength $A_{i}^{(k)}$ in an AC power grid system with an
arbitrary topology (23) decays as a power-law of $\omega$ with an exponent
depending on the graph-theoretic distance $d$ between node $k$ and $i$, i.e.

$A_{i}^{(k)}\overset{\omega\rightarrow\infty}{\sim}\left|\Phi^{[d]}_{ki}\right|\omega^{-2d-1},$
(46)

where $\left|\Phi^{[d]}_{ki}\right|$ is a distance- and node-specific
prefactor but independent on the perturbation frequency.

###### Proof 4.12

To determine the asymptotic behaviour of the response strength
$A_{i}^{(k)}=\omega|R_{i}^{(k)}|$, we first reduce $\mbox{Re}[R_{i}^{(k)}]$
and $\mbox{Im}[R_{i}^{(k)}]$ from (42) and (43) to a common denominator
$M(\omega)$ and obtain the respective numerators $N_{\text{Re}}(\omega)$ and
$N_{\text{Im}}(\omega)$ as polynomials of $\omega$,

$\displaystyle M(\omega)$
$\displaystyle:=\displaystyle\prod_{\ell=0}^{N-1}\left[\left(-\omega^{2}+\lambda^{[\ell]}\right)^{2}+\alpha^{2}\omega^{2}\right],$
(47) $\displaystyle N^{\text{Re}}_{ki}(\omega)$
$\displaystyle:=\displaystyle\sum_{\ell=0}^{N-1}v_{k}^{[\ell]}v_{i}^{[\ell]}\left(-\omega^{2}+\lambda^{[\ell]}\right)Q^{[\ell]}(\omega),\text{
and}$ (48) $\displaystyle N^{\text{Im}}_{ki}(\omega)$
$\displaystyle:=\displaystyle\sum_{\ell=0}^{N-1}v_{k}^{[\ell]}v_{i}^{[\ell]}\left(-\alpha\omega\right)Q^{[\ell]}(\omega)\quad\text{with}$
(49) $\displaystyle Q^{[\ell]}(\omega)$
$\displaystyle:=\prod_{\begin{subarray}{c}\ell^{\prime}=0,\ell^{\prime}\neq\ell\end{subarray}}^{N-1}\left[\left(-\omega^{2}+\lambda^{[\ell^{\prime}]}\right)^{2}+\alpha^{2}\omega^{2}\right].$
(50)

The asymptotic behaviour of $M(\omega)$, $N^{\text{Re}}_{ki}(\omega)$ and
$N^{\text{Im}}_{ki}(\omega)$ as $\omega\rightarrow\infty$ is dominated by the
respective leading terms with the highest power of $\omega$. The denominator
scales asymptotically as

$M(\omega)\overset{\omega\rightarrow\infty}{\sim}\omega^{4N}.$ (51)

For the numerators the leading term depends on a common product
$Q^{[\ell]}(\omega)$. As shown in Appendix A, $Q^{[\ell]}(\omega)$ explicitly
depends on $\lambda^{[\ell]}$ in the form of $Q^{[\ell]}(\omega)$ can be
written in terms of $\lambda^{[\ell]}$ and other variables that are dependent
on the underlying matrix $\mathcal{L}$ but independent of the choice of
$\ell$. Thus, we define

$Q(\lambda^{[\ell]},\omega):=\sum_{j=0}^{2N-2}C^{[j]}(\lambda^{[\ell]})\omega^{4N-4-2j}$
(52)

where the coefficients $C^{[j]}(\lambda^{[\ell]})$ are polynomials in
$\lambda^{[\ell]}$ of degree $j$. Inserting the expression of
$Q(\lambda^{[\ell]},\omega)$ (52) to the numerators (48) and (49), we obtain

$\displaystyle N^{\text{Re}}_{ki}$
$\displaystyle=\sum_{\ell=0}^{N-1}v_{k}^{[\ell]}v_{i}^{[\ell]}\sum_{j=0}^{2N-1}F^{[j]}(\lambda^{[\ell]})\omega^{4N-2-2j}\text{
and }N^{\text{Im}}_{ki}$
$\displaystyle=\sum_{\ell=0}^{N-1}v_{k}^{[\ell]}v_{i}^{[\ell]}\sum_{j=0}^{2N-2}G^{[j]}(\lambda^{[\ell]})\omega^{4N-3-2j}$
(53)

where $F^{[j]}(\lambda^{[\ell]})$ and $G^{[j]}(\lambda^{[\ell]})$ are also
polynomials in $\lambda^{[\ell]}$ of degree $j$ and can be written in terms of
$C^{[j]}(\lambda^{[\ell]})$ as

$\displaystyle F^{[j]}(\lambda^{[\ell]})$
$\displaystyle=\begin{cases}-C^{[j]}(\lambda^{[\ell]})&j=0\\\
-C^{[j]}(\lambda^{[\ell]})+\lambda^{[\ell]}C^{j-1}(\lambda^{[\ell]})&1\leq
j\leq 2N-2\\\ \lambda^{[\ell]}C^{j-1}(\lambda^{[\ell]})&j=2N-1\end{cases}$
(54) $\displaystyle G^{[j]}(\lambda^{[\ell]})$ $\displaystyle=-\alpha
C^{[j]}(\lambda^{[\ell]}),\quad j\in\\{0,1,\cdots,2N-2\\}.$ (55)

Considering the numerators $N^{\text{Re}}_{ki}(\omega)$ and
$N^{\text{Im}}_{ki}(\omega)$ as the $ki$-th elements of numerator matrices
$\mathcal{N}^{\text{Re}}$ and $\mathcal{N}^{\text{Im}}$, we can conveniently
write (53) in a matrix form

$\mathcal{N}^{\text{Re}}=\sum_{j=0}^{2N-1}{\Phi}^{[j]}\omega^{4N-2-2j},\quad\mathcal{N}^{\text{Im}}=\sum_{j=0}^{2N-2}{\Gamma}^{[j]}\omega^{4N-3-2j}$
(56)

with the coefficient matrices
${\Phi}^{[j]}:=\mathcal{V}\mathcal{F}^{[j]}\mathcal{V}^{\text{T}}$ and
${\Gamma}^{[j]}:=\mathcal{V}\mathcal{G}^{[j]}\mathcal{V}^{\text{T}}$. Here
$\mathcal{V}:=(v^{[0]},\cdots,v^{[N-1]})$ and $\mathcal{F}^{[j]}$,
$\mathcal{G}^{[j]}$ are diagonal matrices with
$\mathcal{F}^{[j]}_{ii}:=F^{[j]}(\lambda^{[i]})$ and
$\mathcal{G}^{[j]}_{ii}:=G^{[j]}(\lambda^{[i]})$, respectively. By spelling
out the polynomials

$F^{[j]}(\lambda^{[\ell]})=\sum_{m=0}^{j}f^{[j]}_{m}\cdot(\lambda^{[\ell]})^{m},\quad
G^{[j]}(\lambda^{[i]})=\sum_{m=0}^{j}g^{[j]}_{m}\cdot(\lambda^{[\ell]})^{m}$
(57)

with coefficients $f_{m}^{[j]},g_{m}^{[j]}\in\mathbb{R}$, we can see the
diagonal matrices $\mathcal{F}^{[j]}$ and $\mathcal{G}^{[j]}$ are polynomials
of a diagonal matrix $\Lambda$ with $\Lambda_{ii}:=\lambda^{[i]}$

$\mathcal{F}^{[j]}=\sum_{m=0}^{j}f^{[j]}_{m}\Lambda^{m},\quad\mathcal{G}^{[j]}=\sum_{m=0}^{j}g^{[j]}_{m}\Lambda^{m},$
(58)

so that the coefficient matrices $\Phi^{[j]}$ and $\Gamma^{[j]}$ can be
written as

$\displaystyle{\Phi}^{[j]}$
$\displaystyle=\mathcal{V}\left(\sum_{m=0}^{j}f^{[j]}_{m}\Lambda^{m}\right)\mathcal{V}^{\text{T}}=\sum_{m=0}^{j}f^{[j]}_{m}\left(\mathcal{V}\Lambda^{m}\mathcal{V}^{\text{T}}\right)=\sum_{m=0}^{j}f^{[j]}_{m}\mathcal{L}^{m},$
(59) $\displaystyle{\Gamma}^{[j]}$
$\displaystyle=\mathcal{V}\left(\sum_{m=0}^{j}g^{[j]}_{m}\Lambda^{m}\right)\mathcal{V}^{\text{T}}=\sum_{m=0}^{j}g^{[j]}_{m}\left(\mathcal{V}\Lambda^{m}\mathcal{V}^{\text{T}}\right)=\sum_{m=0}^{j}g^{[j]}_{m}\mathcal{L}^{m},$
(60)

indicating that they are polynomials of the weighted graph Laplacian matrix
$\mathcal{L}$ of degree $j$, as
$\mathcal{V}\Lambda\mathcal{V}^{\text{T}}\equiv\mathcal{L}$. In short, the
numerators $\mathcal{N}^{\text{Re}}$ and $\mathcal{N}^{\text{Im}}_{ki}$ are in
fact polynomials in $\omega$, with its coefficient being also a polynomial in
$\mathcal{L}$,

$\mathcal{N}^{\text{Re}}=\sum_{j=0}^{2N-1}\left(\sum_{m=0}^{j}f^{[j]}_{m}\mathcal{L}^{m}\right)\omega^{4N-2-2j},\quad\mathcal{N}^{\text{Im}}=\sum_{j=0}^{2N-2}\left(\sum_{m=0}^{j}g^{[j]}_{m}\mathcal{L}^{m}\right)\omega^{4N-3-2j}.$
(61)

We note that as the index $j$ increases, the powers of $\omega$ decrease and
the degrees of the polynomials $\Phi^{[j]}(\mathcal{L})$ and
$\Gamma^{[j]}(\mathcal{L})$, i.e. the coefficients of $\omega$, increase. For
sufficiently large perturbation frequency $\omega$, the leading terms in the
numerators which dominate the asymptotic behaviors would be the ones with the
highest degrees of $\omega$ with nonzero coefficients. We know from graph
theory that $(\mathcal{L}^{m})_{ij}\neq 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
$\mathcal{L}$, satisfying $d(i,j)\leq m$. Therefore, the first terms in
$N^{\text{Re}}_{ki}(\omega)$ and $N^{\text{Im}}_{ki}(\omega)$ with the highest
degrees of $\omega$ have exactly zero coefficients, i.e. $\Phi^{[j]}=0$ and
$\Gamma^{[j]}=0$, for all $j<d(k,i)$. Consequently, the leading term in the
numerators are

$N^{\text{Re}}_{ki}(\omega)\overset{\omega\rightarrow\infty}{\sim}\Phi^{[d]}_{ki}\omega^{4N-2-2d},\quad
N^{\text{Im}}_{ki}(\omega)\overset{\omega\rightarrow\infty}{\sim}\Gamma^{[d]}_{ki}\omega^{4N-3-2d}$
(62)

with $d:=d(k,i)$. Together with the asymptotic behavior of the denominator
(51), we obtain

$A_{i}^{(k)}=\omega\left|R_{i}^{(k)}\right|\overset{\omega\rightarrow\infty}{\sim}\omega\cdot\left|\dfrac{\Phi^{[d]}_{ki}\omega^{4N-2-2d}}{\omega^{4N}}\right|=\left|\Phi^{[d]}_{ki}\right|\omega^{-2d-1}.$
(63)

######  Remark 4.13 (Localized response patterns in the high-frequency limit)

4.11 implies that the response amplitude of the grid frequency $A_{i}^{(k)}$
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 $\omega$ for fixed $d$. In the high-frequency limit, a
network’s response to the perturbation is restricted to the perturbed node,
i.e.

$\lim_{\omega\rightarrow\infty}\dfrac{A_{i}^{(k)}}{A_{i}^{(k)}}\overset{\omega\rightarrow\infty}{\sim}\lim_{\omega\rightarrow\infty}\dfrac{\left|\Phi^{[d]}_{ki}\right|\omega^{-2d-1}}{\left|\Phi^{[0]}_{ki}\right|\omega^{-1}}=\lim_{\omega\rightarrow\infty}\left|\Phi^{[d]}_{ki}\right|\omega^{-2d}=\delta_{ki}$
(64)

with $\delta_{ki}$ being the Kronecker delta function.

Figure 2: Topological localization of network responses. For frequencies
larger than all eigenfrequencies and across network types (row 1) the response
amplitudes (45) decays exponentially with shortest-path distance $d$ (row 2)
and algebraically with driving frequency $\omega$ (row 3) (cf. 4.11). Dashed
vertical lines in row 3 indicate the displayed frequency responses in row 2.
Columns display graphs and responses for a random tree (column a), the
topology of the British high voltage transmission grid (column b) and a random
power grid network topology generated according to Schultz et al. (2014)
(column c). Network settings:
$(N,N_{g},P_{g},P_{c},K_{g},K_{c},\alpha)=(264,24,10\text{ s}^{-2},-1\text{
s}^{-2},200\text{ s}^{-2},20\text{ s}^{-2},1\text{ s}^{-1})$ for column a,
$(120,30,39\text{ s}^{-2},-13\text{ s}^{-2},390\text{ s}^{-2},390\text{
s}^{-2},1\text{ s}^{-1})$ for column b, and $(80,20,39\text{ s}^{-2},-13\text{
s}^{-2},390\text{ s}^{-2},390\text{ s}^{-2},1\text{ s}^{-1})$ for column c
333Here in network settings $N_{g}$ is the number of power generating units
with power injection $P_{g}>0$ and the rest units are power consuming units
with $P_{c}<0$. The transmission lines connected to generators have capacity
$K_{g}$ and the rest lines have capacity $K_{c}$..

######  Remark 4.14 (Localized response patterns in networks of multi-
dimensional dynamical systems)

We consider networks of $N$ diffusively coupled identical units governed by
$n$-dimensional dynamics. Each unit $i$ has $n$ state variables
$(x_{i}^{[0]},x_{i}^{[1]}\cdots,x_{i}^{[n-1]})$ and is governed by dynamics

$\displaystyle\begin{cases}\dot{x}_{i}^{[0]}&=x_{i}^{[1]}\\\
\dot{x}_{i}^{[1]}&=x_{i}^{[2]}\\\ \cdots\\\
\dot{x}_{i}^{[n-1]}&=f(x_{i}^{[0]},x_{i}^{[1]}\cdots,x_{i}^{[n-1]})+\displaystyle\sum_{j=1}^{N}g(x_{j}^{[0]}-x_{i}^{[0]}),\end{cases}$

where $f:\mathbb{R}^{N}\rightarrow\mathbb{R}$ and
$g:\mathbb{R}\rightarrow\mathbb{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
$\omega\rightarrow\infty$, we conjecture that the amplitude
$\tilde{A}_{i,m}^{(k)}$ of the sinusoidal response in state variable
$D_{t}^{m}x_{i}$ of unit $i$ is given by

$\tilde{A}_{i,m}^{(k)}\overset{\omega\rightarrow\infty}{\sim}\left|\Psi^{[d]}_{ki}\right|\omega^{-n(d+1)+m},$
(65)

where $|\Psi^{[d]}_{ki}|$ is a distance- and node-specific prefactor but
independent on the perturbation frequency.

The response of a sinusoidally driven damped harmonic oscillator with
$\tilde{A}_{i,m}^{(k)}\overset{\omega\rightarrow\infty}{\sim}\omega^{-2}$ can
be seen as a special case of (65) with $n=2,d=0$ and $m=0$. For networks of
Kuramoto phase oscillators with $n=1$, (65) 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$, (65) reduces to 4.11.

######  Remark 4.15 (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 $\Phi_{ki}^{[d]}$ in (46), is network- and
topology-dependent by definition.

###  4.2 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 Taylor-expand the linear response (28) at node $i$ to a sinusoidal
perturbation at node $k$ in powers of $t$ as

$\Theta_{i}^{(k)}(t)=\sum_{n=0}^{\infty}\dfrac{D_{t}^{n}{\Theta}_{i}^{(k)}(0)}{n!}t^{n}$
(66)

around $t=0$, which is characterized by the time derivatives of the linear
response at $t=0$. Here $D^{n}_{t}:=\frac{\mathrm{d}^{n}}{\mathrm{d}t^{n}}$ is
Euler’s notation for differential operator. The $n$-th order time derivative
of the linear response at $t=0$ is

$D^{n}_{t}\Theta_{i}^{(k)}(0)=\sum_{\ell=0}^{N-1}\dfrac{v_{k}^{[\ell]}v_{i}^{[\ell]}\varepsilon
e^{\imath\varphi}}{-\omega^{2}+\imath\alpha\omega+\lambda^{[\ell]}}\left[\dfrac{\left(\Delta^{[\ell]}_{+}\right)^{n}\left(\Delta^{[\ell]}_{-}-\imath\omega\right)-\left(\Delta^{[\ell]}_{-}\right)^{n}\left(\Delta^{[\ell]}_{+}-\imath\omega\right)}{2\eta^{[\ell]}}+(\imath\omega)^{n}\right].$
(67)

The transient response of a power grid network can thus be estimated by the
first non-zero term in the power series of $t$ (66). 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.

######  Proposition 4.16 (Leading-term approximation of transient response)

The transient response at node $i$ in an AC power grid network (23) to a
sinusoidal perturbation of frequency $\omega$ 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 (28) around $t=0$,

$\Theta_{i}^{(k)}(t)=\dfrac{\varepsilon
e^{\imath\varphi}(-1)^{d}\left(\mathcal{L}^{d}\right)_{ki}}{\left(2d+2\right)!}t^{2d+2}+O(t^{2d+3}).$
(68)

Here $d:=d(k,i)$ denotes the graph-theoretic distance between node $k$ and
node $i$.

###### Proof 4.17

To find out the leading term in the Taylor series of the linear response (66),
we study the summand of the $\ell$-th eigenmode in the derivative of the
linear response $D^{n}_{t}\Theta_{i}^{(k)}(0)$ (67), which we denote as
$v_{k}^{[\ell]}v_{i}^{[\ell]}\varepsilon
e^{\imath\varphi}F_{n}(\lambda^{[\ell]})$ for convenience444Please not that
the functions $F_{n}$ here and $F^{[j]}$ in (54) are different from each other
though both are polynomials in $\lambda^{[\ell]}$.. In the summand, the
function $F_{n}(\lambda^{[\ell]})$ appears to be a division of two polynomials
of $\lambda^{[\ell]}$. In fact, it can be shown that the leading term of
$F_{n}(\lambda^{[\ell]})$ (denoted as
$\mathrm{LT}\left[F_{n}(\lambda^{[\ell]})\right]$), i.e. the term with the
highest order of $\lambda^{[\ell]}$ is

$\displaystyle\mathrm{LT}\left[F_{n}(\lambda^{[\ell]})\right]=\left\\{\begin{array}[]{ll}(-1)^{\frac{n-1}{2}}\left(-\imath\omega+\dfrac{n-1}{2}\alpha\right)\left(\lambda^{[\ell]}\right)^{\frac{n-3}{2}}&\quad\text{if
}n\text{ is odd,}\\\
(-1)^{\frac{n-2}{2}}\left(\lambda^{[\ell]}\right)^{\frac{n-2}{2}}&\quad\text{if
}n\text{ is even.}\end{array}\right.$ (71)

A derivation of the result (71) is given in Appendix B. We note that, for
$n=0$ and $n=1$, $F_{n}(\lambda^{[\ell]})=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 $\Theta_{i}^{(k)}$ and the
frequency response $\dot{\Theta}_{i}^{(k)}$ are zero at the onset of
perturbation. For $n\geq 2$, (71) indicates a monotonic relation between the
degree of $F_{n}(\lambda^{[\ell]})$ as a polynomial of $\lambda^{[\ell]}$ and
the order of derivative $n$: $n=2\deg[F_{n}(\lambda^{[\ell]})]+2$.

As we have shown in Sec. 4.1.3, sums of the form of
$\sum_{\ell=0}^{N-1}v_{k}^{[\ell]}v_{i}^{[\ell]}P^{j}(\lambda^{[\ell]})$,
where $P^{j}(\lambda^{[\ell]})$ represents a general polynomial of
$\lambda^{[\ell]}$ of degree $j$, can be seen as $[P^{j}(\mathcal{L})]_{ki}$,
the $ki$-th element of the matrix $P^{j}(\mathcal{L})$, a polynomial of
$\mathcal{L}$ with degree $j$. Applying this result to the derivatives of the
linear response (67), we find that $D^{n}_{t}\Theta_{i}^{(k)}(0)$ can be
considered as the $ki$-th element of matrix $F_{n}(\mathcal{L})$, a polynomial
of $\mathcal{L}$. The leading term of $D^{n}_{t}\Theta_{i}^{(k)}(0)$ is thus
given by

$\mathrm{LT}\left[D^{n}_{t}\Theta_{i}^{(k)}(0)\right]=\varepsilon
e^{\imath\varphi}\mathrm{LT}\left[F_{n}(\mathcal{L})\right]_{ki},$ (72)

which contains $(\mathcal{L}^{m})_{ki}$ with $m=\frac{n-3}{2}$ if $n$ is odd
and $m=\frac{n-2}{2}$ if $n$ is even. We notice that for a given node pair
$(k,i)$ at distance $d$, we have $(\mathcal{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 (66) 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,

$D^{2d+2}_{t}\Theta_{i}^{(k)}(0)=\varepsilon
e^{\imath\varphi}\mathrm{LT}\left[F_{2d+2}(\mathcal{L})\right]_{ki}=\varepsilon
e^{\imath\varphi}(-1)^{d}(\mathcal{L}^{d})_{ki},$ (73)

because all other terms contain $(\mathcal{L}^{m})_{ki}$ with $m<d$ and thus
vanish. Taken together, the transient linear response near $t=0$ can be
approximated as

$\Theta_{i}^{(k)}(t)=\sum_{n=2d+2}^{\infty}\dfrac{D_{t}^{n}\Theta_{i}^{(k)}(0)}{n!}t^{n}=\dfrac{\varepsilon
e^{\imath\varphi}(-1)^{d}\left(\mathcal{L}^{d}\right)_{ki}}{\left(2d+2\right)!}t^{2d+2}+O(t^{2d+3}).$
(74)

######  Remark 4.18 (Topological factor in perturbation spreading)

The leading-term approximation of the linear response (4.16) 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 $(\mathcal{L}^{d})_{ki}$ in
the leading-term approximation. It satisfies

$\left(\mathcal{L}^{d}\right)_{ki}=\sum_{\mathcal{P}_{k\rightarrow i}^{d}}\
\prod_{(u,v)\in\mathcal{P}_{k\rightarrow i}^{d}}\mathcal{L}_{uv},$ (75)

suggesting that it can be interpreted as the product of the edge weights along
a shortest path $\mathcal{P}_{k\rightarrow i}^{d}$ 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.19 (Scaling behaviours in transient responses)

The leading-term approximation of the transient response (68) reveals two
scaling behaviours as $t\rightarrow 0$ (cf. Fig. 3). First, the transient
response grows algebraically in time with a distance-dependent exponent:
$\Theta_{i}^{(k)}(t)\sim C_{d}t^{2d+2}$. Here $C_{d}:={\varepsilon
e^{\imath\varphi}(-1)^{d}\left(\mathcal{L}^{d}\right)_{ki}}/{\left(2d+2\right)!}$
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\rightarrow 0$.

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
$\Theta_{i}^{(k)}(t)\sim C_{d}t^{2d+2}$ as $t\rightarrow 0$, with time
independent constant $C_{d}$ that depends on signal magnitude, topology, base
operating state and inter-node distance $d$, see (68). 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},\alpha)=(6,3,1\text{ s}^{-2},-1\text{
s}^{-2},10\text{ s}^{-2},10\text{ s}^{-2},1\text{ s}^{-1})$. For the
perturbation signal $(\varepsilon,\omega/2\pi,\varphi)=(1,1\text{ Hz},0\text{
rad})$.

######  Remark 4.20 (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.16) does not depend on the specific choice of network topology.
We underline that the evaluation of the topological factor in perturbation
spreading (Remark 4.18) 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.

###  4.3 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\rightarrow\infty$) and in the
transient stage ($t\rightarrow 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
$I_{\text{res}}:=[\omega_{\text{eigen}}^{[1]},\omega_{\text{eigen}}^{[N-1]}]$,
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_{\text{res}}$ is a task not only of practical significance regarding the
operational safety of power grid systems, but also with a high theoretical
complexity.

######  Definition 4.21 (Dynamic vulnerability index (DVI) for random network
resonances)

In an AC power grid system (23) perturbed by a random fluctuation at node $k$,
characterized by a power spectral density $S(\omega)$ with frequency
components $\omega\in
I_{\mathrm{res}}=[\omega_{\mathrm{eigen}}^{[1]},\omega_{\mathrm{eigen}}^{[N-1]}]$,
the nodal Dynamic Vulnerability Index (DVI) is defined as

$\mathrm{DVI}_{i}^{(k)}:=\int_{I_{\mathrm{res}}}{S(\omega)}^{\frac{1}{2}}\left|\sum_{\ell=0}^{N-1}\dfrac{\imath\omega
v^{[\ell]}_{k}v^{[\ell]}_{i}}{-\omega^{2}+\imath\alpha\omega+\lambda^{[\ell]}}\right|\mathrm{d}\omega.$
(76)

######  Proposition 4.22 (DVI estimates ranking of the nodal all-time-high
steady frequency responses)

Let an AC power grid system (23) be perturbed by a time-dependent fluctuation
at node $k$. A random signal time series is characterized by a power spectral
density $S(\omega)$, i.e., the strength $\varepsilon$ of its frequency
component $\varepsilon e^{\imath(\omega t+\varphi)}$ is frequency dependent
and follows $\varepsilon(\omega)\propto S(\omega)^{\frac{1}{2}}$, and the
corresponding phase $\varphi$ is randomly drawn from the uniform distribution
on $[0,2\pi)$, independently for each realization of the fluctuation time
series. Suppose the fluctuation signal is composed of frequency components
with $\omega\in
I_{\mathrm{res}}=[\omega_{\mathrm{eigen}}^{[1]},\omega_{\mathrm{eigen}}^{[N-1]}]$,
the ranking $\sigma_{\mathrm{ATH}}$ of the nodal all-time-high steady
frequency response magnitude $\max_{t\in[0,T]}|\dot{\Theta}^{(k)}_{i}(t)|$ in
an observation window $T$ approaches the ranking $\sigma_{\mathrm{DVI}}$ of
the DVI defined in Definition 4.21 as the observation time window goes to
infinity, i.e.

$\lim_{T\rightarrow\infty}\sigma_{\mathrm{ATH}}(i)=\sigma_{\mathrm{DVI}}(i)\quad\text{for
all }i\in\\{1,\cdots,N\\}.$ (77)

###### Proof 4.23

To analyse the steady network response, i.e. the response in a time window
$T\rightarrow\infty$, 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 (38). The steady-state response $\dot{\Theta}_{i}(t)$
of the frequency at node $i$ to a perturbation signal $\varepsilon
e^{\imath(\omega t+\varphi)}$ at node $k$ is given by

$\dot{\Theta}_{i}^{(k)}(t)\overset{t\rightarrow\infty}{\sim}\sum_{\ell=0}^{N-1}\dfrac{\imath\varepsilon\omega
v_{k}^{[\ell]}{v}_{i}^{[\ell]}}{-\omega^{2}+\imath\alpha\omega+\lambda^{[\ell]}}e^{\imath(\omega
t+\varphi)}=\imath\varepsilon\omega R_{i}^{(k)}e^{\imath(\omega t+\varphi)},$
(78)

which can be seen as a driven oscillation with a complex amplitude
$\imath\varepsilon\omega R_{i}^{(k)}$. Here $R_{i}^{(k)}$ is the response
factor defined in (37) characterizing the response strength at individual
nodes in a network. The complex amplitude gives rise to an amplitude shift
$|\imath\varepsilon\omega R_{i}^{(k)}|$ and a phase shift
$\arg(\imath\varepsilon\omega R_{i}^{(k)})$, which both are specific to the
perturbation strength $\varepsilon$ and the perturbation frequency $\omega$.

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
$\dot{\Theta}_{i}^{(k)}(\varepsilon,\omega,t)$ (78) over all frequency
components (see also Remark 2.6). Particularly, for perturbation signals which
are characterized by a specific power spectral density (PSD) $S(w)$, the
strength of its frequency component $\varepsilon e^{\imath(\omega t+\varphi)}$
can be expressed in terms of the frequency as $\varepsilon(\omega)\propto
S(\omega)^{\frac{1}{2}}$ while no assumptions is made on the choice of its
phase $\varphi$. 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

$\dot{\Theta}_{i}^{(k)}(S(\omega),t)\overset{t\rightarrow\infty}{\sim}\displaystyle\int_{I_{\text{res}}}\imath
cS(\omega)^{\frac{1}{2}}\omega R_{i}^{(k)}e^{\imath(\omega
t+\varphi)}\mathrm{d}\omega.$ (79)

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 (79) is reached only when the involved frequencies are
finite in number and non-resonant to each other555In this specific context we
adopt the definition of resonance between frequencies in ergodic theory. The
frequencies, as elements of a vector $\bm{\omega}$, are called _resonant_ to
each other if there exists a nonzero integer vector
$\mathbf{m}\cdot\bm{\omega}=0$. Otherwise the frequencies are called non-
resonant to each other. so that all oscillating responses for each frequency
component $\omega$ (78) would eventually align. The largest possible magnitude
of frequency response time series is given by the sum over $\omega$ of the
magnitude of each frequency response $|\imath cS(\omega)^{\frac{1}{2}}\omega
R_{i}^{(k)}|$. For finite time series with $M$ data points and time step
$\Delta T$, its frequency components $\omega_{n}=\frac{2\pi n}{M\Delta T}$
with $n\in\\{1,\cdots,\frac{M}{2}\\}$ 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 $\Delta\omega=2\pi/T$ and thus the more frequency
components exist within the given interval of interest $I_{\text{res}}$. As
$T$ approaches infinity, the order of the resonant frequencies
$\\{\omega_{n}\\}$ becomes sufficiently high so that the alignment of the
phases can be attained at a finite rate dumas_ergodization_1991. Therefore, as
$T\rightarrow\infty$, the all-time-high frequency response approaches the sum
of the frequency-specific response magnitudes, which is proportional to the
DVI given in (76):

$\lim_{T\rightarrow\infty}\max_{t\in[0,T]}\left|\dot{\Theta}^{(k)}_{i}(S(\omega),t)\right|=\int_{I_{\text{res}}}cS(\omega)^{\frac{1}{2}}\left|\imath\omega
R_{i}^{(k)}\right|\mathrm{d}\omega=c\,\text{DVI}_{i}^{(k)}.$ (80)

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 (80) does not play a role in the ranking and we
finally arrive at

$\lim_{T\rightarrow\infty}\sigma_{\mathrm{ATH}}(i)=\sigma_{\mathrm{DVI}}(i)\quad\text{for
all }i\in\\{1,\cdots,N\\}.$ (81)

######  Remark 4.24 (Generalization of DVI and convergence time)

The DVI defined in (76) 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_{\text{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 80 and 81) 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.

##  5 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.

###  5.1 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\rightarrow\infty$ (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 $\alpha$, 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 (28). 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 $\mathcal{L}$. Through
asymptotic analyses, either with respect to time $t$ or the perturbation
frequency $\omega$, 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 $\omega$ approaches its corresponding limit (see 4.11 and
4.16). As the contribution of each eigenmode contains the “overlap factor”
$v_{i}^{[\ell]}v_{k}^{[\ell]}$ of the perturbed node $k$ and the responding
node $i$, the powers of the Laplacian eigenvalue $(\lambda^{[\ell]})^{m}$ that
is involved in the dominating term in each eigenmode $\ell$ translates to
elements of the power of the Laplacian matrix $(\mathcal{L}^{m})_{ki}$ through
the summation over all $N$ eigenmodes $\ell\in\\{0,\cdots,N-1\\}$.
Furthermore, with the help of the result from graph theory that
$(\mathcal{L}^{m})_{ki}\equiv 0$ if $m\in\mathbb{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 $\omega$ (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.18 to estimate the contribution of
higher order terms $O(2d+3)$ in the Taylor series (68) 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)\sim t^{d}$ where $d$ is the shortest-path distance
between perturbed node and the node $i$ the response in measured at.

###  5.2 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 (28) 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 $\dot{\bm{\Theta}}$ (78) to a single-frequency sinusoidal
oscillation $D_{k}(t)$ at a given node $k$ is a mapping
$f_{G,\bm{\theta}^{*}}:\mathbb{R}\rightarrow\mathbb{R}^{N}$ with
$f_{G,\bm{\theta}^{*}}$ depending on the underlying network topology $G$ and
network’s base state $\bm{\theta}^{*}$ 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 $\varepsilon$, but becomes
further dependent of the frequency $\omega$ through the PSD $S(\omega)$ of the
signal, i.e. $\varepsilon(\omega)=S(\omega)^{\frac{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_{\text{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.

###  5.3 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 (27 and 28) 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 $\bm{\theta}^{*}$ (29), the deviation
of the nonlinear coupling terms $\sin(\theta_{j}-\theta_{i})$ for all edges
$(i,j)\in E$ to their values $\sin(\theta^{*}_{j}-\theta^{*}_{i})$ at the
fixed point are represented by the first-order approximations
$\cos(\theta^{*}_{j}-\theta^{*}_{i})(\Theta_{j}-\Theta_{i})$, vectorized as
the term $-\mathcal{L}\bm{\Theta}$ in (29). For power grid systems working at
a stable operation state without any transmission line overloaded
($|\theta_{j}^{*}-\theta_{i}^{*}|\leq\frac{\pi}{2}$ for all edges $(i,j)\in
E$, see 2.2), the linear approximation breaks down only when the system is
driven far enough from the fixed point $\bm{\theta}^{*}$ and goes close to the
point where one of the lines $(i,j)$ is fully loaded, i.e.
$\sin(\theta_{j}-\theta_{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].

##  6 Conclusions and outlook

###  6.1 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 time-dependent 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.

###  6.2 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
$\frac{\mathrm{d}g_{ij}}{\mathrm{d}(x_{j}-x_{i})}$ 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 (9). 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 $\mathcal{J}$ in (3) is diagonal or almost diagonal (such as in the
Jordan normal forms of $\mathcal{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 (37)
or even (28) 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.

###### Acknowledgements.

The authors acknowledge the support from the Deutsche Forschungsgemeinschaft
(DFG; German Research Foundation) under Germany’s Excellence Strategy
EXC-2068-390729961- Cluster of Excellence Physics of Life, the Center for
Advancing Electronics at TU Dresden, by the German Federal Ministry for
Research and Education (BMBF grants no. 03SF0472F and no. 03EK3055F), the
National Natural Science Foundation of China (grant no. 12161141016), Shanghai
Municipal Science and Technology Major Project (grant no. 2021SHZDZX0100), and
the Shanghai Municipal Commission of Science and Technology Project (grant
nos. 18ZR1442000 and 19511132101).

##  Appendix A Proof of Equation 52

###### Proposition A.1

Given $\omega>0,\alpha>0$ and $0=\lambda^{[0]}<\cdots<\lambda^{[N-1]}$ as
defined in subsection 3.2, the product $Q^{[\ell]}(\omega)$

$Q^{[\ell]}(\omega):=\prod_{\begin{subarray}{c}\ell^{\prime}=0,\ell^{\prime}\neq\ell\end{subarray}}^{N-1}\left[\left(-\omega^{2}+\lambda^{[\ell^{\prime}]}\right)^{2}+\alpha^{2}\omega^{2}\right]$
(82)

that appears in the numerators of the real part and of the imaginary part of
the nodal response strength (40) explicitly depends on $\lambda^{[\ell]}$ and
can be expressed as

$Q(\lambda^{[\ell]},\omega)=\sum_{j=0}^{2N-2}C^{[j]}(\lambda^{[\ell]})\omega^{4N-4-2j}$
(83)

with the coefficient $C^{[j]}(\lambda^{[\ell]})$ is a polynomial of
$\lambda^{[\ell]}$ with degree $j$.

###### Proof A.2

To prove the proposition, we first rewrite the factors in $Q^{[\ell]}(w)$ by
ordering the terms according to the degree of $\omega$:

$Q^{[\ell]}(w)=\prod_{\begin{subarray}{c}\ell^{\prime}=0,\ell^{\prime}\neq\ell\end{subarray}}^{N-1}\left[\omega^{4}+\left(\alpha^{2}-2\lambda^{[\ell^{\prime}]}\right)\omega^{2}+\left(\lambda^{[\ell^{\prime}]}\right)^{2}\right]=:\prod_{\begin{subarray}{c}\ell^{\prime}=0,\ell^{\prime}\neq\ell\end{subarray}}^{N-1}\sum_{m=1}^{3}r_{m}(\lambda^{[\ell^{\prime}]})\omega^{2(3-m)}.$
(84)

According to the distributivity of multiplication over addition, it is clear
from (84) that each term in $Q^{[\ell]}(w)$, a polynomial of $\omega$ with
degree $4N-4$, can be seen as the product of three factors:
$\prod_{\ell^{\prime}\in{s_{1}}}r_{1}(\lambda^{[\ell^{\prime}]})\omega^{4}$,
$\prod_{\ell^{\prime}\in{s_{2}}}r_{2}(\lambda^{[\ell^{\prime}]})\omega^{2}$
and $\prod_{\ell^{\prime}\in{s_{3}}}r_{3}(\lambda^{[\ell^{\prime}]})$, where
sets $s_{1}$, $s_{2}$, and $s_{3}$ have $a:=|s_{1}|$, $b:=|s_{2}|$ and
$c:=|s_{3}|$ elements respectively and together form a partition
$P_{\ell}(a,b,c)$ of the set of the indices of the $N-1$ eigenmodes
$S_{\ell}:=\\{0,...,N-1\\}\backslash\\{\ell\\}$. Here $a,b,c\in\mathbb{N}_{0}$
and satisfy $a+b+c=N-1$. Since $r_{1}$, $r_{2}$ and $r_{3}$ are polynomials of
$\lambda^{[\ell^{\prime}]}$ with degree $0$, $1$ and $2$, a term in
$Q^{[\ell]}(w)$ with $\omega^{4a+2b}$ would have a coefficient involving a
multiplication of $2b+c$ Laplacian eigenvalues $\lambda^{[\ell^{\prime}]}$
with $\ell^{\prime}\in S_{\ell}$. Denoting $j=b+2c$, we can write the
coefficient of the term with degree $4a+2b=4N-4-2j$ as

$C^{[\ell]}_{j}=\sum_{\begin{subarray}{c}a+b+c=N-1\\\
b+2c=j\end{subarray}}\sum_{P_{\ell}(a,b,c)}\prod_{p\in
s_{2}}\left(\alpha^{2}-2\lambda^{[p]}\right)\prod_{q\in
s_{3}}\left(\lambda^{[q]}\right)^{2},$ (85)

which is a sum over all possible partitions $P_{\ell}(a,b,c)$ satisfying
$a+b+c=N-1$ and $b+2c=j$.

In the following we show that the coefficient $C^{[\ell]}_{j}$ is a polynomial
of $\lambda^{[\ell]}$ with degree $j$, i.e.
$\deg[C^{[\ell]}_{j}(\lambda^{[\ell]})]=j$. For convenience of notation in the
proof, we define the sum of coefficients involving both $r_{2}$ and $r_{3}$
over $s_{2}\in\binom{S_{\ell}}{b}$ and $s_{3}\in\binom{S_{\ell}\backslash
s_{2}}{c}$, i.e. all possible partitions $P_{\ell}(a,b,c)$ of $S_{\ell}$ as

$Y^{[\ell]}_{b,c}:=\sum_{s_{2}\in\binom{S_{\ell}}{b},s_{3}\in\binom{S_{\ell}\backslash
s_{2}}{c}}\prod_{p\in s_{2}}\left(\alpha^{2}-2\lambda^{[p]}\right)\prod_{q\in
s_{3}}\left(\lambda^{[q]}\right)^{2}.$ (86)

Here $\binom{S_{\ell}}{b}$ denotes all possible $b$-subsets of $S_{\ell}$.
Similarly, we define the sum of the coefficients over all possible partitions
of $S:=\\{0,...,N-1\\}$ as

$Y_{b,c}:=\sum_{s_{2}\in\binom{S}{b},s_{3}\in\binom{S\backslash
s_{2}}{c}}\prod_{p\in s_{2}}\left(\alpha^{2}-2\lambda^{[p]}\right)\prod_{q\in
s_{3}}\left(\lambda^{[q]}\right)^{2}.$ (87)

In case $b=0$ or $c=0$, the corresponding product is omitted. It is clear that
$Y_{b,c}$, including special cases $Y_{b,0}$ and $Y_{0,c}$ are constants
independent of $\lambda^{[\ell]}$. Using definition (86) and (87), we can
write $C_{j}^{[\ell]}$ in 85 as
$\sum_{\begin{subarray}{c}a+b+c=N-1,b+2c=j\end{subarray}}Y^{[\ell]}_{b,c}$. To
prove A.1, we only need to prove $Y^{[\ell]}_{b,c}$ is a polynomial of
$\lambda^{[\ell]}$ with degree $j=b+2c$, i.e.

$\deg\left[Y^{[\ell]}_{b,c}(\lambda^{[\ell]})\right]=b+2c.$ (88)

Now we show (88) in three steps. All subproofs are given by mathematical
induction.

  1. [ (P1)]

  2. Step 1:

First, we show that the sum of $r_{2}$-related factors over
$s_{2}\in\binom{S_{\ell}}{b}$ is a polynomial of $\lambda^{[\ell]}$ with
degree $b$. That is, $\deg[Y^{[\ell]}_{(b,0)}(\lambda^{[\ell]})]=b$.

    1. [(P1)]

    2. (a)

For $b=1$, we have $Y^{[\ell]}_{1,0}=Y_{1,0}-(\alpha^{2}-2\lambda^{[\ell]})$,
which is a polynomial of $\lambda^{[\ell]}$ with degree $1$ since $Y_{1,0}$ is
a constant independent of $\lambda^{[\ell]}$.

    3. (b)

If the statement holds for $b=n-1$, i.e.
$\deg[Y^{[\ell]}_{n-1,0}(\lambda^{[\ell]})]=n-1$, then for $b=n$ we have
$Y^{[\ell]}_{n,0}=Y_{n,0}-\left(\alpha^{2}-2\lambda^{[\ell]}\right)Y^{[\ell]}_{n-1,0}$,
satisfying $\deg[Y^{[\ell]}_{n,0}(\lambda^{[\ell]})]=n$.

  3. Step 2:

Second, we show that the sum of $r_{3}$-related factors over
$s_{3}\in\binom{S_{\ell}}{c}$ is a polynomial of $\lambda^{[\ell]}$ with
degree $2c$. That is, $\deg[Y^{[\ell]}_{(0,c)}(\lambda^{[\ell]})]=2c$.

    1. [(P1)]

    2. (a)

For $c=1$, we have $Y^{[\ell]}_{0,1}=Y_{0,1}-(\lambda^{[\ell]})^{2}$, which is
a polynomial of $\lambda^{[\ell]}$ with degree $2$ since $Y_{0,1}$ is a
constant independent of $\lambda^{[\ell]}$.

    3. (b)

If the statement holds for $c=n-1$, i.e.
$\deg[Y^{[\ell]}_{0,n-1}(\lambda^{[\ell]})]=2n-2$, then for $c=n$ we have
$Y^{[\ell]}_{0,n}=Y_{0,n}-\left(\lambda^{[\ell]}\right)^{2}Y^{[\ell]}_{0,n-1}$,
satisfying $\deg[Y^{[\ell]}_{0,n}(\lambda^{[\ell]})]=2n$.

  4. Step 3:

Finally, we show that the sum of coefficients involving both $r_{2}$ and
$r_{3}$ over $s_{2}\in\binom{S_{\ell}}{b}$ and
$s_{3}\in\binom{S_{\ell}\backslash\\{s_{2}\\}}{c}$, i.e. all possible
partitions $P_{\ell}(a,b,c)$, is a polynomial of $\lambda^{[\ell]}$ with
degree $j=b+2c$. That is, equation (88).

    1. [(P1)]

    2. (a)

For $b=1,c=1$, we have
$Y^{[\ell]}_{1,1}=Y_{1,1}-(\alpha^{2}-2\lambda^{[\ell]})Y^{[\ell]}_{0,1}-(\lambda^{[\ell]})^{2}Y^{[\ell]}_{1,0}$,
which is a polynomial of $\lambda^{[\ell]}$ with degree $3$ since
$\deg[Y_{(1,1)}(\lambda^{[\ell]})]=0$,
$\deg[Y^{[\ell]}_{(0,1)}(\lambda^{[\ell]})]=2$ and
$\deg[Y^{[\ell]}_{(1,0)}(\lambda^{[\ell]})]=1$.

    3. (b)

If the statement holds for $b=m-1,c=n-1$, i.e.
$\deg[Y^{[\ell]}_{m-1,n-1}(\lambda^{[\ell]})]=m+2n-3$, then for $b=m,c=n$ we
have

$\displaystyle Y^{[\ell]}_{m,n}=$ $\displaystyle
Y_{m,n}-\left(\alpha^{2}-2\lambda^{[\ell]}\right)Y^{[\ell]}_{m-1,n}-\left(\lambda^{[\ell]}\right)^{2}Y^{[\ell]}_{m,n-1}$
$\displaystyle=$ $\displaystyle
Y_{m,n}-\left(\alpha^{2}-2\lambda^{[\ell]}\right)\left(Y_{m-1,n}-\left(\lambda^{[\ell]}\right)^{2}Y^{[\ell]}_{m-1,n-1}\right)$
$\displaystyle-\left(\lambda^{[\ell]}\right)^{2}\left(Y_{m,n-1}-\left(\alpha^{2}-2\lambda^{[\ell]}\right)Y^{[\ell]}_{m-1,n-1}\right).$

Taking into account that $Y_{m,n}$, $Y_{m-1,n}$ and $Y_{m,n-1}$ all have
degree $0$, we can easily see that
$\deg[Y^{[\ell]}_{m,n}(\lambda^{[\ell]})]=m+2n$, meaning the statement also
holds for $b=m$ and $c=n$.

##  Appendix B Proof of Equation 71

###### Proposition B.1

The function

$F_{n}(\lambda^{[\ell]}):=\dfrac{1}{-\omega^{2}+\imath\alpha\omega+\lambda^{[\ell]}}\left[\dfrac{\left(\Delta^{[\ell]}_{+}\right)^{n}\left(\Delta^{[\ell]}_{-}-\imath\omega\right)-\left(\Delta^{[\ell]}_{-}\right)^{n}\left(\Delta^{[\ell]}_{+}-\imath\omega\right)}{2\eta^{[\ell]}}+(\imath\omega)^{n}\right]$
(89)

that appears in the $n$-th order derivative of the linear response at $t=0$
(67) has a leading term with respect to $\lambda^{[\ell]}$

$\displaystyle\mathrm{LT}\left[F_{n}(\lambda^{[\ell]})\right]=\left\\{\begin{array}[]{ll}(-1)^{\frac{n-1}{2}}\left(-\imath\omega+\frac{n-1}{2}\alpha\right)\left(\lambda^{[\ell]}\right)^{\frac{n-3}{2}}&\quad\text{if
}n\text{ is odd,}\\\
\left(-\lambda^{[\ell]}\right)^{\frac{n-2}{2}}&\quad\text{if }n\text{ is
even.}\end{array}\right.$ (92)

Here $\Delta^{[\ell]}_{\pm}:=-\alpha/2\pm\eta^{[\ell]}$,
$\eta^{[\ell]}:=\sqrt{\alpha^{2}/4-\lambda^{[\ell]}}$ with $\alpha>0$,
$0=\lambda^{[0]}<\cdots<\lambda^{[N-1]}$ and $\omega>0$, $n\in\mathbb{N}$,
$n\geq 2$.

###### Proof B.2

Using the relation $\Delta^{[\ell]}_{+}\Delta^{[\ell]}_{-}=\lambda^{[\ell]}$
we rewrite the function under study as

$F_{n}(\lambda^{[\ell]})=\dfrac{\lambda^{[\ell]}f_{n-1}(\lambda^{[\ell]})-\imath\omega
f_{n}(\lambda^{[\ell]})+(\imath\omega)^{n}}{-\omega^{2}+\imath\alpha\omega+\lambda^{[\ell]}}$
(93)

with

$f_{n}(\lambda^{[\ell]}):=\dfrac{1}{2\eta^{[\ell]}}\left[(\Delta^{[\ell]}_{+})^{n}-(\Delta^{[\ell]}_{-})^{n}\right].$
(94)

It is clear from (93) that the leading term of $F_{n}(\lambda^{[\ell]})$
depends on the leading term of $f_{n}(\lambda^{[\ell]})$ as

$\displaystyle\mathrm{LT}\left[F_{n}(\lambda^{[\ell]})\right]$
$\displaystyle=\frac{\mathrm{LT}\left[\lambda^{[\ell]}f_{n-1}(\lambda^{[\ell]})-\imath\omega
f_{n}(\lambda^{[\ell]})+(\imath\omega)^{n}\right]}{\mathrm{LT}\left[-\omega^{2}+\imath\alpha\omega+\lambda^{[\ell]}\right]}$
(95)
$\displaystyle=\frac{1}{\lambda^{[\ell]}}\mathrm{LT}\left[\lambda^{[\ell]}f_{n-1}(\lambda^{[\ell]})-\imath\omega
f_{n}(\lambda^{[\ell]})\right].$ (96)

Please note that $(\Delta^{[\ell]}_{+})^{n}$ and $(\Delta^{[\ell]}_{-})^{n}$
in $f_{n}(\lambda^{[\ell]})$ are a complex conjugate pair, since
$\eta^{[\ell]}$ is imaginary under the low dissipation of power grid systems
(see Remark 3.7). Therefore we have
$(\Delta^{[\ell]}_{-})^{n}=\overline{(\Delta^{[\ell]}_{+})^{n}}$ which leads
to
$f_{n}(\lambda^{[\ell]})=\text{Im}\,(\Delta^{[\ell]}_{+})^{n}/\sqrt{\lambda^{[\ell]}-\alpha^{2}/4}$.
Now proving (92) boils down to determining the leading term of
$(\Delta^{[\ell]}_{+})^{n}$. In the following we use mathematical induction to
show that leading term of the real part and the imaginary part of
$(\Delta^{[\ell]}_{+})^{n}$ follows

$\displaystyle\mathrm{LT}\left[\text{Re}\,(\Delta^{[\ell]}_{+})^{n}\right]$
$\displaystyle=\left\\{\begin{array}[]{ll}(-1)^{\frac{n+1}{2}}\frac{n}{2}\alpha\left(\lambda^{[\ell]}\right)^{\frac{n-1}{2}}&\quad\text{if
}n\text{ is odd,}\\\
\left(-\lambda^{[\ell]}\right)^{\frac{n}{2}}&\quad\text{if }n\text{ is
even;}\\\ \end{array}\right.$ (99)
$\displaystyle\mathrm{LT}\left[f_{n}(\lambda^{[\ell]})\right]=\mathrm{LT}\left[\dfrac{\text{Im}\,(\Delta^{[\ell]}_{+})^{n}}{\sqrt{\lambda^{[\ell]}-\alpha^{2}/4}}\right]$
$\displaystyle=\left\\{\begin{array}[]{ll}\left(-\lambda^{[\ell]}\right)^{\frac{n-1}{2}}&\quad\text{if
}n\text{ is odd,}\\\
(-1)^{\frac{n}{2}}\frac{n}{2}\alpha\left(\lambda^{[\ell]}\right)^{\frac{n-2}{2}}&\quad\text{if
}n\text{ is even.}\\\ \end{array}\right.$ (102)

  1. (a)

For $n=2$ and $n=3$, we can easily verify (99) and (102) by spelling out
$(\Delta^{[\ell]}_{+})^{n}$:

$\displaystyle\mathrm{LT}\left[\text{Re}\,(\Delta^{[\ell]}_{+})^{2}\right]=\mathrm{LT}\left[-\lambda^{[\ell]}+\tfrac{1}{2}\alpha^{2}\right]=-\lambda^{[\ell]},$
$\displaystyle\mathrm{LT}\left[\text{Re}\,(\Delta^{[\ell]}_{+})^{3}\right]=\mathrm{LT}\left[\tfrac{3}{2}\alpha\lambda^{[\ell]}-\tfrac{1}{4}\alpha^{2}-\tfrac{1}{4}\alpha^{3}\right]=\tfrac{3}{2}\alpha\lambda^{[\ell]},$
$\displaystyle\mathrm{LT}\left[f_{2}(\lambda^{[\ell]})\right]=\mathrm{LT}\left[-\alpha\right]=-\alpha,\quad\mathrm{LT}\left[f_{3}(\lambda^{[\ell]})\right]=\mathrm{LT}\left[-\lambda^{[\ell]}+\alpha^{2}\right]=-\lambda^{[\ell]}.$

  2. (b)

Now we show that (99) and (102) hold for $n+1$ if they hold for $n$, no matter
$n$ is odd or even. The leading term of
$\mbox{Re}\,(\Delta^{[\ell]}_{+})^{n+1}$ and $f_{n+1}(\lambda^{[\ell]})$ can
be expressed in terms of the leading term of
$\mbox{Re}\,(\Delta^{[\ell]}_{+})^{n}$ and $f_{n}(\lambda^{[\ell]})$ as
following

$\displaystyle\mathrm{LT}\left[\mbox{Re}\,(\Delta^{[\ell]}_{+})^{n+1}\right]$
$\displaystyle=\mathrm{LT}\left[(-\tfrac{1}{2}\alpha)\mathrm{LT}\left[\mbox{Re}\,(\Delta^{[\ell]}_{+})^{n}\right]-\mathrm{LT}\left[f_{n}(\lambda^{[\ell]})\right](\lambda^{[\ell]}-\tfrac{1}{4}\alpha^{2})\right],$
(103) $\displaystyle\mathrm{LT}\left[f_{n+1}(\lambda^{[\ell]})\right]$
$\displaystyle=\mathrm{LT}\left[\mathrm{LT}\left[\mbox{Re}\,(\Delta^{[\ell]}_{+})^{n}\right]+(-\tfrac{1}{2}\alpha)\mathrm{LT}\left[f_{n}(\lambda^{[\ell]})\right]\right].$
(104)

In case $n$ is odd, we have

$\displaystyle\mathrm{LT}\left[\mbox{Re}\,(\Delta^{[\ell]}_{+})^{n+1}\right]$
$\displaystyle=\mathrm{LT}\left[(-1)^{\frac{n+3}{2}}\tfrac{n}{4}\alpha^{2}(\lambda^{[\ell]})^{\frac{n-1}{2}}+(-\lambda^{[\ell]})^{\frac{n-1}{2}}(\lambda^{[\ell]}-\tfrac{1}{4}\alpha^{2})\right]$
$\displaystyle=(-\lambda^{[\ell]})^{\frac{n+1}{2}},\text{ and}$ (105)
$\displaystyle\mathrm{LT}\left[f_{n+1}(\lambda^{[\ell]})\right]$
$\displaystyle=\mathrm{LT}\left[(-1)^{\frac{n+1}{2}}\tfrac{n}{2}\alpha(\lambda^{[\ell]})^{\frac{n-1}{2}}-\tfrac{1}{2}\alpha(-\lambda^{[\ell]})^{\frac{n-1}{2}}\right]$
$\displaystyle=(-1)^{\frac{n+1}{2}}\tfrac{n+1}{2}\alpha(\lambda^{[\ell]})^{\frac{n-1}{2}}$
(106)

In case $n$ is even, we have

$\displaystyle\mathrm{LT}\left[\mbox{Re}\,(\Delta^{[\ell]}_{+})^{n+1}\right]$
$\displaystyle=\mathrm{LT}\left[(-\lambda^{[\ell]})^{\frac{n}{2}}(-\tfrac{1}{2}\alpha)+(-1)^{\frac{n}{2}}\dfrac{n}{2}\alpha(\lambda^{[\ell]})^{\frac{n-2}{2}}(\lambda^{[\ell]}-\tfrac{1}{4}\alpha^{2})\right]$
$\displaystyle=(-1)^{\frac{n+2}{2}}\tfrac{n+1}{2}\alpha(\lambda^{[\ell]})^{\frac{n}{2}},\text{
and}$ (107) $\displaystyle\mathrm{LT}\left[f_{n+1}(\lambda^{[\ell]})\right]$
$\displaystyle=\mathrm{LT}\left[(-\lambda^{[\ell]})^{\frac{n}{2}}+(-1)^{\frac{n+2}{2}}\tfrac{n+1}{2}\alpha(\lambda^{[\ell]})^{\frac{n-2}{2}}\right]=(-\lambda^{[\ell]})^{\frac{n}{2}}.$
(108)

The results (105, 106, 107, 108) agree with the statement (99) and (102).

Combining results (96) and (102), we arrive at the leading term of
$F_{n}(\lambda^{[\ell]})$ as

$\displaystyle\mathrm{LT}\left[F_{n}(\lambda^{[\ell]})\right]$
$\displaystyle=\dfrac{1}{\lambda^{[\ell]})}\mathrm{LT}\left[(-1)^{\frac{n-1}{2}}\tfrac{n-1}{2}\alpha(\lambda^{[\ell]})^{\frac{n-1}{2}}-\imath\omega(-\lambda^{[\ell]})^{\frac{n-1}{2}}\right]$
$\displaystyle=(-1)^{\frac{n-1}{2}}\left(-\imath\omega+\tfrac{n-1}{2}\alpha\right)(\lambda^{[\ell]})^{\frac{n-3}{2}}\quad\text{if
}n\text{ is odd, and}$ (109)
$\displaystyle\mathrm{LT}\left[F_{n}(\lambda^{[\ell]})\right]$
$\displaystyle=\dfrac{1}{\lambda^{[\ell]})}\mathrm{LT}\left[(-1)^{\frac{n-2}{2}}(\lambda^{[\ell]})^{\frac{n}{2}}-\imath\omega(-1)^{\frac{n}{2}}\tfrac{n}{2}\alpha(\lambda^{[\ell]})^{\frac{n-2}{2}}\right]$
$\displaystyle=(-\lambda^{[\ell]})^{\frac{n-2}{2}}\quad\text{if }n\text{ is
even.}$ (110)

## References

  * (1)
  * Acebrón et al. (2005) Acebrón, J. A., Bonilla, L. L., Pérez Vicente, C. J., Ritort, F. and Spigler, R. (2005).  The Kuramoto model: A simple paradigm for synchronization phenomena, Reviews of Modern Physics 77(1): 137–185.
  * Anvari et al. (2016) Anvari, M., Lohmann, G., Wächter, M., Milan, P., Lorenz, E., Heinemann, D., Tabar, M. R. R. and Peinke, J. (2016).  Short term fluctuations of wind and solar power systems, New Journal of Physics 18(6): 063027.
  * Bapat (2010) Bapat, R. B. (2010).  Graphs and matrices, Universitext, Springer ; Hindustan Book Agency, London ; New York : New Delhi.  OCLC: ocn455828013.
  * Cammi and Mennucci (1999) Cammi, R. and Mennucci, B. (1999).  Linear response theory for the polarizable continuum model, The Journal of Chemical Physics 110(20): 9877–9886.
  * Casagrande (2006) Casagrande, V. (2006).  Synchronization, waves, and turbulence in systems of interacting chemical oscillators.  Publisher: Technische Universität Berlin.
  * Coletta et al. (2018) Coletta, T., Bamieh, B. and Jacquod, P. (2018).  Transient Performance of Electric Power Networks Under Colored Noise, 2018 IEEE Conference on Decision and Control (CDC), IEEE, Miami Beach, FL, pp. 6163–6167.
  * Dörfler et al. (2013) Dörfler, F., Chertkov, M. and Bullo, F. (2013).  Synchronization in complex oscillator networks and smart grids, Proceedings of the National Academy of Sciences 110(6): 2005–2010.
  * ENTSO-E (2021) ENTSO-E (2021).  Continental Europe Synchronous Area Separation on 8 January 2021: Interim Report, Technical report, European Network of Transmission System Operators for Electricity.  https://eepublicdownloads.azureedge.net/clean-documents/Publications/Position%20papers%20and%20reports/entso-e_CESysSep_interim_report_210225.pdf.
  * Filatrella et al. (2008) Filatrella, G., Nielsen, A. H. and Pedersen, N. F. (2008).  Analysis of a power grid using a Kuramoto-like model, The European Physical Journal B 61(4): 485–491.