# Supplementary Material to:
Resonant Solitary States in Complex Networks

Jakob Niehues  Potsdam Institute for Climate Impact Research, Telegrafenberg
A56, 14473 Potsdam, Germany  Humboldt-Universität zu Berlin, Department of
Physics, Newtonstraße 15, 12489 Berlin, Germany  Technische Universität
Berlin, ER 3-2, Hardenbergstrasse 36a, 10623 Berlin, Germany  S. Yanchuk
Potsdam Institute for Climate Impact Research, Telegrafenberg A56, 14473
Potsdam, Germany  Humboldt-Universität zu Berlin, Department of Physics,
Newtonstraße 15, 12489 Berlin, Germany  University College Cork, School of
Mathematical Sciences, Western Road, Cork, T12 XF62, Ireland  R. Berner
Humboldt-Universität zu Berlin, Department of Physics, Newtonstraße 15, 12489
Berlin, Germany  J. Kurths  Potsdam Institute for Climate Impact Research,
Telegrafenberg A56, 14473 Potsdam, Germany  Humboldt-Universiät zu Berlin,
Department of Physics, Newtonstraße 15, 12489 Berlin, Germany  F. Hellmann
Potsdam Institute for Climate Impact Research, Telegrafenberg A56, 14473
Potsdam, Germany  M. Anvari  Potsdam Institute for Climate Impact Research,
Telegrafenberg A56, 14473 Potsdam, Germany  Fraunhofer Institute for
Algorithms and Scientific Computing, 53757 Sankt Augustin, Germany

This supplementary material contains a detailed derivation of the central
self-consistent equation in the main article. Moreover, we provide details on
the examples studied.

###### Contents

  1. I Derivation of the self-consistent equation
    1. I.1 Preliminaries
      1. I.1.1 Network modes
      2. I.1.2 Eigenvector localization
      3. I.1.3 Partition of the system
    2. I.2 Averaging approach
    3. I.3 Linear response
    4. I.4 Validity
    5. I.5 Solitary nodes with several neighbors
    6. I.6 Heterogeneous system parameters
    7. I.7 Several solitary nodes
  2. II Evaluation of the self-consistent equation
    1. II.1 The minimal effective model
    2. II.2 The star graph
    3. II.3 Complex synthetic power grids

##  I Derivation of the self-consistent equation

###  I.1 Preliminaries

We start with a network of $N$ Kuramoto oscillators with inertia (swing
equations),

$m_{i}\ddot{\varphi}_{i}=P_{i}-\alpha_{i}\dot{\varphi}_{i}-\sum_{j=1}^{N}\kappa_{ij}\sin(\varphi_{i}-\varphi_{j}).$
(S1)

A necessary condition for the existence of a steady state is that the driving
powers are balanced, that is, $\sum_{i=1}^{N}P_{i}=0$. As long as the sum of
damping coefficients is not zero, this can always be achieved by shifting to a
rotating frame of reference with
$\omega_{\text{ref}}=\frac{\sum_{i=1}^{N}P_{i}}{\sum_{i=1}^{N}\alpha_{i}}$.
Without loss of generality, we will assume $\sum_{i=1}^{N}P_{i}=0$ from here
on. Alternatively, this can be written as a set of $2N$ first-order equations,

$\displaystyle\dot{\varphi}_{i}$ $\displaystyle=$ $\displaystyle\omega_{i},$
(S2) $\displaystyle m_{i}\dot{\omega}_{i}$ $\displaystyle=$ $\displaystyle
P_{i}-\alpha_{i}\omega_{i}-\sum_{j=1}^{N}\kappa_{ij}\sin(\varphi_{i}-\varphi_{j}).$
(S3)

The steady state $\varphi_{i}^{*}$, i.e., $\dot{\varphi}_{i}=0$, is given by
the power flow equations,

$P_{i}=\sum_{j=1}^{N}\kappa_{ij}\sin(\varphi_{i}^{*}-\varphi_{j}^{*}).$ (S4)

####  I.1.1 Network modes

The steady state induces a weighted steady-state graph Laplacian

$\mathcal{L}_{ij}=\begin{cases}\begin{aligned} \sum_{j,j\neq
i}^{N}\kappa_{ij}\cos\left(\varphi_{i}^{*}-\varphi_{j}^{*}\right)&\quad\text{iff
$i=j$},\\\
-\kappa_{ij}\cos\left(\varphi_{i}^{*}-\varphi_{j}^{*}\right)&\quad\text{iff
$i\neq j$.}\end{aligned}\end{cases}$ (S5)

It appears in the linearized dynamics with Jacobian $\mathcal{J}$,

$\frac{{\text{d}}}{{\text{d}}t}\begin{bmatrix}\mathds{1}&0\\\
0&\mathcal{M}\end{bmatrix}\begin{bmatrix}\delta\varphi\\\
\delta\omega\end{bmatrix}=\mathcal{J}\begin{bmatrix}\delta\varphi\\\
\delta\omega\end{bmatrix}=\begin{bmatrix}0&\mathds{1}\\\
-\mathcal{L}&-\mathcal{A}\end{bmatrix}\begin{bmatrix}\delta\varphi\\\
\delta\omega\end{bmatrix}.$ (S6)

where $\delta\varphi$ is a vector with entries $\varphi_{i}-\varphi_{i}^{*}$,
analogously $\delta\omega$ is a vector with entries
$\omega_{i}-\omega_{i}^{*}=\omega_{i}$, and
$\mathcal{M}\coloneqq\text{diag}(m_{i})$, and
$\mathcal{A}\coloneqq\text{diag}(\alpha_{i})$. $\mathcal{L}$ is positive semi-
definite, see reference [7]. If the graph is connected, there is exactly one
zero eigenvalue. The other eigenvalues are real and positive, and can be
ordered
$0=\Lambda^{\left[1\right]}<\Lambda^{\left[2\right]}<\ldots<\Lambda^{\left[N\right]}$.
For $\mathcal{M}=\mathds{1}$ and $\mathcal{A}=\alpha\mathds{1}$ 111Note that
for homogeneous $m_{i}$ and $\alpha_{i}$, one can always be rescaled to 1 by
rescaling time., the eigenvalues of $\mathcal{J}$ are given by

$\sigma_{\pm}^{\left[\ell\right]}=-\frac{\alpha}{2}\pm\sqrt{\frac{\alpha^{2}}{4}-\Lambda^{\left[\ell\right]}},$
(S7)

see reference [26]. Since their real parts are negative by assumption 222We
only treat networks with linearly stable steady states., except for one zero
eigenvalue, the synchronous state is linearly stable as small perturbations
are damped exponentially with time. The imaginary parts of
$\sigma_{\pm}^{\left[\ell\right]}$ are the eigenmodes of the network, or
network modes, $\pm\omega^{\left[\ell\right]}$, that describe oscillations
around the fixed point. They are given by $\omega^{[1]}=0$ and

$\omega^{\left[\ell\right]}=\sqrt{\Lambda^{\left[\ell\right]}-\frac{\alpha^{2}}{4}}.$
(S8)

for $\ell=2,\ldots,N$. We assume that
$\Lambda^{\left[\ell\right]}\gg\alpha^{2}/4$ for $\ell\geq\ell_{\text{th}}$,
where $\ell_{\text{th}}\geq 2$ is some threshold. Therefore,
$\omega^{[\ell]}\approx\sqrt{\Lambda^{[\ell]}}$ for
$\ell\geq\ell_{\text{th}}$.

The components $V_{i}^{[\ell]}$ of the corresponding eigenvectors
$\vec{V}^{[\ell]}$ give us the nodal excitation of $\varphi_{i}$ of network
mode $\ell$ at node location $i$. The zero eigenmode $\omega^{[1]}=0$
corresponds to the uniform eigenvector $V_{i}^{[1]}=1/\sqrt{N}\quad\forall i$.
It represents the symmetry of Eq. (S1) under a constant global phase shift. It
does not describe an oscillation, but the uniform movement of all phases.
Therefore, it is also called bulk mode.

####  I.1.2 Eigenvector localization

Here, we provide an estimate for the eigenvalue of a highly concentrated
steady-state Laplacian eigenvector. It is a typical feature of the networks at
hand and will become important in the derivation of the self-consistent
equation. Let us assume that a network mode is localized at node $i$, or as
one could say, its Laplacian eigenvector $\vec{V}^{[\ell]}$ has a strong
support there. We quantify this by

$\left|\frac{V_{j}^{[\ell]}}{V_{i}^{[\ell]}}\right|<\varepsilon_{i}\quad\forall
j\in\mathcal{N}_{i},$ (S9)

where $\varepsilon_{i}$ is a small parameter, and we denote the neighbors of
node $i$ with $\mathcal{N}_{i}$. When evaluated explicitly for component $i$,
the eigenvector equation yields

$\Lambda^{[\ell]}V_{i}^{[\ell]}=\mathcal{L}_{ii}V_{i}^{[\ell]}+\sum_{j\in\mathcal{N}_{i}}\mathcal{L}_{ij}V_{j}^{[\ell]}=\mathcal{L}_{ii}V_{i}^{[\ell]}-\sum_{j\in\mathcal{N}_{i}}\kappa_{ij}\cos\left(\varphi_{i}^{*}-\varphi_{j}^{*}\right)V_{j}^{[\ell]}.$
(S10)

This gives us tight bounds for the value of $\Lambda^{[\ell]}$,

$\left|\Lambda^{[\ell]}-\mathcal{L}_{ii}\right|=\left|\sum_{j\in\mathcal{N}_{i}}\kappa_{ij}\cos\left(\varphi_{i}^{*}-\varphi_{j}^{*}\right)\frac{V_{j}^{[\ell]}}{V_{i}^{[\ell]}}\right|\leq\mathcal{L}_{ii}\varepsilon_{i}.$
(S11)

For $\cos\left(\varphi_{i}^{*}-\varphi_{j}^{*}\right)\approx 1$, this reduces
to

$\left|\Lambda^{[\ell]}-K_{i}^{*}d_{i}\right|\leq
K_{i}^{*}d_{i}\varepsilon_{i},$ (S12)

where $K_{i}^{*}$ is the average coupling strength of the links
$\\{(i,j)|j\in\mathcal{N}_{i}\\}$ with coupling strenghts $\\{\kappa_{ij}\\}$,
and $d_{i}=|\mathcal{N}_{i}|$ is the degree of node $i$. We conclude that a
highly localized network mode comes with a Laplacian eigenvalue that is
closely connected to the degree of the node where the network mode is
localized.

####  I.1.3 Partition of the system

For the analysis of solitary states, we view the original system as two
interacting subsystems, namely, the solitary oscillator $z$ and the rest,
which we denote $S=\\{1,\ldots,N\\}\setminus z$, and which forms a
synchronized cluster. We assume that node $z$ is of degree one with neighbor
$k$. As we show below, the case of several neighbors can be treated
analogously due to the additive coupling.

Node $z$ has a different mean frequency than the synchronized cluster, i.e.,

$\varphi_{z}(t)=\omega_{s}t+\varepsilon\psi_{z}(t),$ (S13)

where $\omega_{s}$ is the time-averaged frequency of node $z$,
$\omega_{s}\coloneqq\langle\dot{\varphi}_{z}\rangle$, $\varepsilon$ a small
parameter, and $\langle\psi_{z}(t)\rangle=0$ by assumption. We assume
homogeneous inertia $m$ and damping constants $\alpha$ for the synchronized
cluster $S$. We denote the energy flow from node $i$ to node $j$ as

$p_{ij}\coloneqq\kappa_{ij}\sin(\varphi_{i}-\varphi_{j})=-p_{ji}.$ (S14)

For now, we assume that the solitary oscillator has only one neighbor, $k$,
without loss of generality. As we show later on, the results can be
generalized to several neighbors in a straightforward way due to the additive
nature of the coupling via the links from $z$ to $k\in\mathcal{N}_{z}$.

For the solitary oscillator, we have

$m_{z}\ddot{\varphi}_{z}=P_{z}-\alpha_{z}\dot{\varphi}_{z}-p_{zk}.$ (S15)

We assume that the solitary state is sufficiently close to being periodic with
a period that is a multiple of $T_{s}\coloneqq 2\pi/\omega_{s}$, hence the
energy flows are oscillating functions in time. Regarding the coupling of the
two subsystems, we separate the mean energy flow and the oscillating, mean-
zero energy flow,

$p_{zk}(t)\coloneqq\langle p_{zk}\rangle+p_{zk,\text{osc}}(t).$ (S16)

With the ansatz from Eq. (S13), time averaging of Eq. (S15) yields a self-
consistent equation for the solitary mean frequency,

$0=Z(\omega_{s})\coloneqq P_{z}-\alpha_{z}\omega_{s}-\langle p_{zk}\rangle.$
(S17)

In order to evaluate Eq. (S17), we need an approximation for $\langle
p_{zk}\rangle$ as a function of $\omega_{s}$, thus closing the equation. By
assuming $\omega_{s}$ to be generally slowly varying, we can provide the
following heuristic dynamical interpretation. The self-consistency function
$Z(\omega_{s})$ is the average change in $m_{z}\dot{\varphi}_{z}$ over a
period, which needs to be zero for a stable, periodic solitary state. In other
words, solutions $\omega_{s}$ to $Z(\omega_{s})=0$ are fixed points
$\omega_{s}=v^{*}=\text{const.}$ of the dynamical system $\dot{v}=Z(v)$ with
state variable $v$. Hence, the sign of $\mathrm{d}Z/\mathrm{d}\omega_{s}$
gives us the linear stability of the zeros of $Z(\omega_{s})$, providing us
with a a heuristic for the stability of the corresponding periodic orbit with
$\langle\dot{\varphi}_{z}\langle=\omega_{s}$.

To determine $\langle p_{zk}\rangle$ as a function of $\omega_{s}$, system
parameters and topology, we first need to find an explicit approximation for
the trajectory of $\varphi_{k}(t)$, as Eq. (S1) can not be solved
analytically.

We start by considering $\langle p_{zk}\rangle$ and $p_{zk,\text{osc}}(t)$ one
at a time. First, we need to account for the power-imbalance in $S$. Eq. (S17)
shows that non-zero $\omega_{s}$ causes a deviation in $p_{zk}$ from its value
in the synchronous state, $P_{z}$. Furthermore, Fig. 1 in the main article
shows that the mean coupling between the solitary node and $S$ is maximal in
the synchronous state, intermediate for intermediate solitary frequency, and
approximately zero for natural solitary states. Hence, in a solitary state,
where the power in the synchronized cluster is unbalanced, the oscillators in
$S$ perform a slow collective counter-rotation to the solitary node, which can
be observed numerically. Adding up Eq. (S1) for the synchronized cluster
yields

$m\frac{\mathrm{d}}{\mathrm{d}t}\sum_{i\in
S}\dot{\varphi}_{i}=-\alpha\sum_{i\in S}\dot{\varphi}_{i}-P_{z}+p_{zk}.$ (S18)

We start by considering only the mean energy flow $\langle p_{zk}\rangle$
between node $z$ and the synchronized cluster. Then, the power in the
synchronized cluster is unbalanced, and the solution of Eq. (S18)
exponentially converges towards

$\sum_{i\in S}\dot{\varphi}_{i}=\frac{-P_{z}+\langle p_{zk}\rangle}{\alpha}.$
(S19)

Therefore, the synchronization frequency of the synchronized cluster is not
zero as in the original system, but, on average, lies at

$\omega_{\text{sync}}\coloneqq\frac{-P_{z}+\langle
p_{zk}\rangle}{\alpha(N-1)}=-\frac{\alpha_{z}}{\alpha}\frac{\omega_{s}}{N-1},$
(S20)

according to Eq. (S17). Note that for constant ratios $m_{i}/\alpha_{i}$ for
all nodes including $z$, the $\alpha_{i}$-weighted sum over all frequencies
exponentially decays towards zero. Therefore, if $m_{z}/\alpha_{z}=m/\alpha$,
we have $\alpha_{z}\omega_{s}=-\alpha(N-1)\omega_{\text{sync}}$.

To be able to linearize the dynamics, we need the proper coordinate frame,
that eliminates the collective rotation of $S$. Therefore, we define the
corotating solitary frequency,
$\omega_{c}\coloneqq\omega_{s}-\omega_{\text{sync}}$. If we shift the phases
in the synchronized cluster by $\omega_{\text{sync}}t$, the synchronous
periodic orbit becomes a fixed point, which induces a steady-state graph
Laplacian $L$. We assume that it has an eigenvector $\vec{v}^{[r]}$ with a
high concentration at node $k$. The corresponding eigenvalue is
$\lambda^{[r]}$. We observe that this is a typical situation for dense root
nodes, which are neighbors of dense sprout nodes. As we show in the
preliminaries, the frequency of a highly localized mode is closely connected
to the degree of the respective node. We remark that, if such a localized
eigenmode does not exist, we can just use the largest eigenmode for the
rescaling in the subsequent derivation. Later on, we see that all modes
contribute to $\langle p_{zk}\rangle$. The magnitude of their contribution
depends on their localization at node $k$ and their frequency.

###  I.2 Averaging approach

Here we derive an explicit expression for Eq. (S17) using the averaging
technique. The first step is to rewrite system (S1) so that the fast and slow
variables are split with a suitable small parameter $\varepsilon$. We consider
the case $\alpha_{i}=\alpha$ and $m_{i}=m$ for all $i\in
S=\\{1,\ldots,N\\}\setminus z$ in the synchronized cluster, and rewrite (S1)
by separating the solitary node $z$ as follows

$\displaystyle m\ddot{\varphi}_{i}$
$\displaystyle=P_{i}-\alpha\dot{\varphi}_{i}-\sum_{j=1}^{N}\kappa_{ij}\sin(\varphi_{i}-\varphi_{j}),\quad
i\in S,$ (S21) $\displaystyle m_{z}\ddot{\varphi}_{z}$
$\displaystyle=P_{z}-\alpha_{z}\dot{\varphi}_{z}-\kappa_{kz}\sin(\varphi_{z}-\varphi_{k}),\quad
k\in\mathcal{N}_{z}.$ (S22)

We first rescale parameters and the time variable333Note that this implies an
inverse rescaling of frequencies. using the following ansatz

$t=\frac{\sqrt{m}}{\nu}t^{\prime},\qquad
P_{i}=\frac{P^{\prime}_{i}\alpha\nu}{\sqrt{m}},\qquad\kappa_{ij}=\frac{\kappa^{\prime}_{ij}\alpha\nu}{\sqrt{m}},$
(S23)

and omitting the prime for the dimensionless quantities (e.g. $t^{\prime}\to
t$) thereafter. Here $\nu$ will be specified later. We also introduce the
small dimensionless parameter

$\varepsilon\coloneqq\frac{\alpha}{\sqrt{m}\nu}\ll 1.$ (S24)

With the rescaling (S23)–(S24), system (S1) transforms into

$\displaystyle\ddot{\varphi}_{i}=$
$\displaystyle\varepsilon\left[P_{i}-\dot{\varphi}_{i}-\sum_{j=1}^{N}\kappa_{ij}\sin(\varphi_{i}-\varphi_{j})\right],\quad
i\in S,$ (S25) $\displaystyle\ddot{\varphi}_{z}=$
$\displaystyle\frac{\varepsilon
m}{m_{z}}\left[P_{z}-\frac{\alpha_{z}}{\alpha}\dot{\varphi}_{z}-\kappa_{zk}\sin(\varphi_{z}-\varphi_{k})\right].$
(S26)

If $\varepsilon$ is small, and the terms multiplied by it on the right-hand
side are of order one, then the averaging theorem tells us that we have
approximately freely rotating oscillators. If we set $\nu=1$, small
$\varepsilon$ means we have a system with very high inertia and low damping.
For such a system, this is indeed a good approximate solution. However, it
does not describe any coupling between the oscillators.

To obtain the solitary and resonant states we want to describe, we want to
consider the case where $\varepsilon$ is small in part due to $\nu$, and $\nu$
is related to the coupling. If we fix $\nu^{2}\varepsilon\approx 1$ then
$\kappa^{\prime}\approx\kappa$, and the right-hand side in (S26) is indeed
small. However, there are $\nu^{2}$ such that the sum in the right-hand side
of (S25) is not small. We will see that these are related to the eigenvalues
of the Laplacian of the synchronized cluster. Thus, we proceed to linearize
the system near the synchronized cluster.

With $\delta_{ij}$ as the Kronecker delta, we define shifted power constants

$\hat{P}_{i}\coloneqq P_{i}-\alpha\omega_{\text{sync}}+\delta_{ik}\langle
p_{zk}\rangle=P_{i}+\frac{\alpha_{z}\omega_{s}}{N-1}+\delta_{ik}\langle
p_{zk}\rangle.$ (S27)

They are balanced, $\sum_{i\in S}\hat{P}_{i}=0$. For the phase angles of the
synchronized cluster, we make the ansatz

$\varphi_{i}(t)=\omega_{\text{sync}}t+\vartheta_{i}^{*}+\vartheta_{i}(t),$
(S28)

where the steady-state of the synchronized cluster with only $\langle
p_{zk}\rangle$ as input is given by $\vartheta_{i}(t)=0$, and

$\hat{P}_{i}=\sum_{j\in
S}\kappa_{ij}\sin\left(\vartheta_{i}^{*}-\vartheta_{j}^{*}\right)$ (S29)

are the power flow equations for the steady-state in the synchronized cluster.
Equation (S28) can be interpreted as a transformation into the average
corotating frame of the synchronized cluster. We define

$\vartheta_{z}(t)\coloneqq\varphi_{z}(t)-\omega_{\text{sync}}t-\vartheta_{k}^{*}=\omega_{c}t+\varepsilon\psi_{z}(t)-\vartheta_{k}^{*}.$
(S30)

The function $\psi_{z}(t)$ is zero with only $\langle p_{zk}\rangle$ as
coupling. The constant phase shift $\vartheta_{k}^{*}$ in $\vartheta_{z}(t)$
can be chosen by applying a translation in time and turns out to be a
convenient choice later.

Next, we consider also the oscillating part of the energy flow between the two
subsystems, $p_{zk,\text{osc}}(t)$. With Eq. (S27), the dynamics Eqs. (S25)
and (S26) for the corotating coordinates read

$\displaystyle\ddot{\vartheta}_{i}$ $\displaystyle=$
$\displaystyle\varepsilon\left[\hat{P}_{i}-\dot{\vartheta}_{i}-\sum_{j=1}^{N}\kappa_{ij}\sin(\vartheta_{i}^{*}+\vartheta_{i}-\vartheta_{j}^{*}-\vartheta_{j})\right],\qquad\text{for
$i\in S\setminus k$,}$ (S31) $\displaystyle\ddot{\vartheta}_{k}$
$\displaystyle=$
$\displaystyle\varepsilon\left[\hat{P}_{k}-\dot{\vartheta}_{k}-\sum_{j\in
S}\kappa_{kj}\sin\left(\vartheta_{k}^{*}+\vartheta_{k}-\vartheta_{j}^{*}-\vartheta_{j}\right)+\kappa_{zk}\sin(\vartheta_{z}-\vartheta_{k})-\langle
p_{zk}\rangle\right],$ (S32) $\displaystyle\ddot{\vartheta}_{z}$
$\displaystyle=$ $\displaystyle\frac{\varepsilon
m}{m_{z}}\left[P_{z}-\frac{\alpha_{z}}{\alpha}(\dot{\vartheta}_{z}+\omega_{\text{sync}})-\kappa_{zk}\sin(\vartheta_{z}-\vartheta_{k})\right].$
(S33)

We remark that due to the definition of $\hat{P}_{k}$, the last two terms in
Eq. (S32) together have zero mean, hence we have successfully eliminated any
constant offset in the perturbation induced by the coupling $p_{zk}$, which
would cause a rotation. We observe that $\vartheta_{k}$ is being harmonically
driven. Due to the additional nonlinear perturbation
$\vartheta_{k}\sin\vartheta_{z}$, we can not simply apply linear response.
Therefore, in the following we apply the averaging theorem to treat this term
and show that it averages out in Eq. (S32). We show later on that one obtains
the same result as with the averaging technique, if one neglects this term
right away, but in Eq. (S32) only, where it is only a small contribution to
the overall coupling of node $k$. It should not be ignored in Eq. (S33), as it
gives a relevant contribution to $\langle p_{zk}\rangle$ in Eq. (S17).

Next, we approximate the full energy flow by linearization with respect to
$\vartheta_{k}$,

$p_{zk}\approx\kappa_{zk}\left[\sin(\vartheta_{z})-\vartheta_{k}\cos\vartheta_{z}\right].$
(S34)

Using Eq. (S29), we linearize Eqs. (S31) - (S33) with respect to
$\vartheta_{i}(t)$ for $i\in S$),

$\displaystyle\ddot{\vartheta}_{i}$ $\displaystyle=$
$\displaystyle\varepsilon\left[-\dot{\vartheta}_{i}-\sum_{j\in
S}\left(\vartheta_{i}-\vartheta_{j}\right)\kappa_{ij}\cos\left(\vartheta_{i}^{*}-\vartheta_{j}^{*}\right)\right]\qquad\text{for
$i\in S\setminus k$,}$ (S35) $\displaystyle\ddot{\vartheta}_{k}$
$\displaystyle=$
$\displaystyle\varepsilon\left[-\dot{\vartheta}_{k}-\sum_{j\in
S}\left(\vartheta_{k}-\vartheta_{j}\right)\kappa_{kj}\cos\left(\vartheta_{k}^{*}-\vartheta_{j}^{*}\right)+\kappa_{zk}\left(\sin\vartheta_{z}-\vartheta_{k}\cos\vartheta_{z}\right)-\langle
p_{zk}\rangle\right],$ (S36) $\displaystyle\ddot{\vartheta}_{z}$
$\displaystyle=$ $\displaystyle\frac{\varepsilon
m}{m_{z}}\left[P_{z}-\frac{\alpha_{z}}{\alpha}(\dot{\vartheta}_{z}+\omega_{\text{sync}})-\left(\kappa_{zk}\sin\vartheta_{z}-\vartheta_{k}\cos\vartheta_{z}\right)\right].$
(S37)

The steady state in the corotating frame of $S$ induces a graph Laplacian for
the synchronized cluster,

$L_{ij}\coloneqq\begin{cases}\begin{aligned} \sum_{j\in S\setminus
i}\kappa_{ij}\cos\left(\vartheta_{i}^{*}-\vartheta_{j}^{*}\right)&\quad\text{iff
$i=j$},\\\
-\kappa_{ij}\cos\left(\vartheta_{i}^{*}-\vartheta_{j}^{*}\right)&\quad\text{iff
$i\neq j$.}\end{aligned}\end{cases}$ (S38)

for $i,j\in S$. We remark that it is reminiscent of the Laplacian of the
original system, $\mathcal{L}$, if one removes all entries associated to node
$z$, i.e., row and column $z$, as well as the contribution to
$\mathcal{L}_{kk}$.

Let $\vec{\vartheta}$ be the vector with entries $\vartheta_{i}$ for $i\in S$.
We can now write Eqs. (S35) and (S36) in vectorized form,

$\displaystyle\ddot{\vec{\vartheta}}$ $\displaystyle=$
$\displaystyle\varepsilon\left(-\dot{\vec{\vartheta}}-L\vec{\vartheta}+\vec{f}\right),$
(S39) $\displaystyle\vec{f}$ $\displaystyle\coloneqq$
$\displaystyle\hat{e}_{k}\left(\kappa_{zk}\sin\omega_{c}t-\kappa_{zk}\vartheta_{k}\cos\omega_{c}t-\langle
p_{zk}\rangle\right),$ (S40)

where $\hat{e}_{k}$ is the Cartesian basis vector in $k$-direction, which has
components $\delta_{ik}$. We have used the zeroth order in $\varepsilon$ for
$\vartheta_{z}$, that is, $\vartheta_{z}\approx\omega_{c}t$.

We assume that there is a network mode with associated eigenvalue
$\lambda^{[r]}$ that is highly localized at node $k$. As shown above, this
implies that the associated frequency is closely connected to $d_{k}$, the
degree of node $k$. We set $\nu^{2}=\lambda^{[r]}$ to measure the systems
frequencies in units of $\sqrt{\lambda^{[r]}}$. Therefore, the above
approximations are most valid for typical root-sprout configurations with
highly connected root nodes.

Since $\kappa^{\prime}_{ij}=\frac{\kappa_{ij}}{\varepsilon\lambda^{[r]}}$ in
Eq. (S23), the Laplacian got rescaled like
$L^{\prime}=\frac{L}{\varepsilon\lambda^{[r]}}$ (renaming $L^{\prime}$ to $L$
thereafter). We now rescale the Laplacian again, such that we measure it in
units of $\lambda^{[r]}$ and make the $\epsilon$-scaling explicit, that is,
$L^{\prime}=\frac{L^{\prime\prime}}{\varepsilon}$, thus substituting
$L=\frac{L^{\prime\prime}}{\varepsilon}$ in Eq. (S39), and renaming
$L^{\prime\prime}$ to $L$ afterward,

$\ddot{\vec{\vartheta}}=\varepsilon\left(-\dot{\vec{\vartheta}}-+\vec{f}\right)-L\vec{\vartheta}.$
(S41)

We diagonalize the system by transforming into the eigenbasis
$\left\\{\vec{v}^{[\ell]}\right\\}$ of $L$,

$\displaystyle\vec{\xi}$ $\displaystyle\coloneqq$ $\displaystyle
Q^{T}\vec{\vartheta},$ (S42) $\displaystyle Q_{ij}$ $\displaystyle\coloneqq$
$\displaystyle v^{[j]}_{i},$ (S43) $\displaystyle Q^{T}LQ$
$\displaystyle\coloneqq$
$\displaystyle\text{diag}(\lambda^{[\ell]})\coloneqq\Gamma,$ (S44)

with the eigenvalues $\left\\{\lambda^{[\ell]}\right\\}$ ordered by magnitude
from $\ell=0$ to $N-1$. The $\xi_{\ell}$ are the excitations of network modes
corresponding to $\omega^{[\ell]}\approx\sqrt{\lambda^{[\ell]}}$ with their
nodal excitations $\vec{\vartheta}$ distributed across the network according
to $\vec{v}^{[\ell]}$.

Eq. (S41) becomes

$\displaystyle\ddot{\vec{\xi}}$ $\displaystyle=$
$\displaystyle\varepsilon\left(-\dot{\vec{\xi}}+\vec{g}\right)-\Gamma\vec{\xi},$
(S45)

where $\vec{g}\coloneqq Q^{T}\vec{f}$. Equation (S45) has the form of the
perturbation of the Hamiltonian system $\ddot{\vec{\xi}}=-\Gamma\vec{\xi}$,
where the perturbation is due to the friction term $\dot{\vec{\xi}}$ and
$\vartheta_{z}$ in $\vec{g}$. The appropriate perturbative approach in this
case is averaging [28], where one should first write the system for the slowly
varying amplitudes of the network modes, i.e., the amplitudes of the periodic
solutions of the unperturbed system. These zeroth order solutions in
$\varepsilon$ of are

$\xi_{\ell}(t)=x_{\ell}\cos\left(\sqrt{\lambda^{[\ell]}}\cdot
t\right)-z_{\ell}\sin\left(\sqrt{\lambda^{[\ell]}}\cdot t\right),$ (S46)

where $\vec{x}$ and $\vec{z}$ are constants specified by initial conditions.
Since the zeroth order system is conservative, they describe a set of periodic
orbits on tori in phase space. For $\ell=1$, we have $\lambda^{[1]}=0$ and an
additional constant plus linear motion in time. However, this uniform motion
of the synchronized cluster is already accounted for in Eq. (S13) via a global
phase shift in the original system, hence we can set it to zero here.

Following reference [28] (theorem 4.1.1, example II), for the first-order
solution, we make the ansatz

$\vec{\xi}=\vec{x}(t)\cos\omega_{c}t-\vec{z}(t)\sin\omega_{c}t,$ (S47)

where $\vec{x}(t)$ and $\vec{z}(t)$ are now time-dependent amplitude
variables. Due to coupling, damping, and resulting resonance between the
solitary rotation and the oscillation of the synchronized cluster, the
$\xi_{\ell}$ slowly move on the solution manifold of the zeroth order system,
and the frequency of their collective oscillation is shifted from the
Laplacian eigenvalue to $\omega_{s}$.

We write the equations of motion, Eqs. (S45), as a system of first-order
differential equations,

$\frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix}\vec{\xi}\\\
\vec{\zeta}\end{pmatrix}\coloneqq\begin{pmatrix}\vec{\zeta}\\\
\varepsilon\left(-\vec{\zeta}+\vec{g}\right)-\Gamma\vec{\xi}\end{pmatrix},$
(S48)

and transform into a corotating frame with the amplitudes as new variables,

$\begin{pmatrix}\vec{x}\\\ \vec{z}\end{pmatrix}\coloneqq
A(t)\begin{pmatrix}\vec{\xi}\\\
\vec{\zeta}\end{pmatrix}\coloneqq\begin{pmatrix}\mathds{1}\cos\omega_{c}t&-\frac{\mathds{1}}{\omega_{c}}\sin\omega_{c}t\\\
-\mathds{1}\sin\omega_{c}t&-\frac{\mathds{1}}{\omega_{c}}\cos\omega_{c}t\end{pmatrix}\begin{pmatrix}\vec{\xi}\\\
\vec{\zeta}\end{pmatrix},$ (S49)

using the invertible van der Pol transformation $A(t)$, see reference [28]
(theorem 4.1.1, example II). The equations of motion for the amplitude
variables are

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix}\vec{x}\\\
\vec{z}\end{pmatrix}$ $\displaystyle=$
$\displaystyle\frac{\mathds{1}\omega_{c}^{2}-\Gamma}{\omega_{c}}\begin{pmatrix}\vec{z}\sin^{2}\omega_{c}t-\vec{x}\cos\omega_{c}t\sin\omega_{c}t\\\
\vec{z}\cos\omega_{c}t\sin\omega_{c}t-\vec{x}\cos^{2}\omega_{c}t\end{pmatrix}$
(S50)
$\displaystyle-\frac{\varepsilon}{\omega_{c}}\begin{pmatrix}\sin\omega_{c}t\left(\omega_{c}\vec{x}\sin\omega_{c}t+\omega_{c}\vec{z}\cos\omega_{c}t+\vec{g}\right)\\\
\cos\omega_{c}t\left(\omega_{c}\vec{x}\sin\omega_{c}t+\omega_{c}\vec{z}\cos\omega_{c}t+\vec{g}\right)\end{pmatrix}.$

To apply the averaging theorem, we need the right-hand side to be of order
$\varepsilon$. For the minimal effective model that we discuss below, the
synchronized cluster only has a single frequency ($\Gamma$ is scalar), and we
can just assume a close resonance, $\omega_{c}^{2}-\Gamma=O(\varepsilon)$. For
the networked case, we need to assume one of the following two possibilities
for all modes. First, if a mode is in resonance, we have
$\omega_{c}^{2}-\lambda=O(\varepsilon)$. Second, if the driving frequency is
not close to the eigenmode’s intrinsic frequency $\omega^{[\ell]}$, say
$\omega_{c}^{2}-\lambda=O(1)$, the mode $z_{\ell}$ and $x_{\ell}$ will be only
weakly excited.

This can be understood heuristically by the observation that in the absence of
a resonant excitation, the response will be of the order of the driving force.
There are two independent reasons for the driving force to be small, the weak
localization of most modes at node $k$, and the fact that the driving comes
with a factor $\varepsilon$.

To make this more explicit, consider the transfer function of Eq. (S45),

$\displaystyle T(s)$ $\displaystyle=$ $\displaystyle(s^{2}+\varepsilon
s+\Gamma)^{-1}$ (S51) $\displaystyle\tilde{\xi}(s)$ $\displaystyle=$
$\displaystyle\varepsilon T(s)\tilde{g}(s),$ (S52)

where $\tilde{\xi}(s)$ and $\tilde{g}(s)$ are the Laplace transforms of
$\vec{\xi}(t)$ and $\vec{g}(t)$, respectively. Absence of resonance implies
that $\tilde{g}(s)$ has no support near $\Gamma$, and thus $T(s)g(s)$ is never
of order $\varepsilon^{-1}$.

We note that this is an additional assumption to example II of the averaging
theorem as in reference [28], which does not treat the networked case $|S|>1$.
As we show later on, this assumption is self-consistent with the solutions it
produces, and with empirically observed behavior of the full system in
simulations.

We define the diagonal matrix $\Delta$ by

$\mathds{1}\omega_{c}^{2}-\Gamma\coloneqq\varepsilon\Delta.$ (S53)

Since the system Eq. (S50) is slow moving it is of the right form to apply the
averaging theorem in reference [28] (theorem 4.1.1 therein). It states that a
system of this form, where the right-hand side is periodic and of order
$\varepsilon$, can be approximated at order $O(\varepsilon)$ by an autonomous
system which is obtained by time-averaging the right-hand side for fixed
(average) $\vec{x}$ and $\vec{z}$. The periodicity is maintained in a
$O(\varepsilon^{2})$ term, which we neglect here. The solutions of both
systems are $\varepsilon$-close for time scales up to $O(\varepsilon^{-1})$
(in dimensionless time). Translating this back by inverting Eq. (S23), this
amounts to time scales of the order $m/\alpha$, which is typically much larger
than $T_{s}=2\pi/\omega_{s}$. Further, hyperbolic fixed points of the averaged
system possess a unique hyperbolic orbit in the original system
$\varepsilon$-close to the fixed point for all $\varepsilon$ up to some
$\varepsilon_{0}$. We make use of the last statement to determine an
approximation for the (average) orbit $\vartheta_{k}$ from the unique stable
fixed point of the averaged system.

By averaging Eq. (S50) over time, we obtain the linear system

$\frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix}\vec{x}\\\
\vec{z}\end{pmatrix}=J\begin{pmatrix}\vec{x}\\\
\vec{z}\end{pmatrix}-\vec{h}\coloneqq-\frac{\varepsilon}{2\omega_{c}}\begin{pmatrix}\mathds{1}\omega_{c}&-\Delta\\\
\Delta&\mathds{1}\omega_{c}\end{pmatrix}\begin{pmatrix}\vec{x}\\\
\vec{z}\end{pmatrix}-\vec{h},$ (S54)

where the components of the inhomogeneity vector are given by

$h_{\ell}=\begin{cases}\begin{aligned}
\frac{\varepsilon\kappa_{zk}}{2\omega_{c}}v_{k}^{[\ell]}&\quad\text{for
$\ell=1,\ldots,N-1$}\\\ 0&\quad\text{for
$\ell=N,\ldots,2N-2$}.\end{aligned}\end{cases}$ (S55)

We note that the nonlinear term in $\vec{g}$ is averaged out, thereby
explaining why we can obtain the same result by neglecting it and applying
linear response, as shown below. The eigenvalues of the Jacobian $J$ are given
by

$j_{\pm}^{[\ell]}=-\frac{\varepsilon}{2}\pm
i\frac{\varepsilon\Delta_{\ell}}{2\omega_{c}},$ (S56)

where $\Delta_{\ell}$ are diagonal elements of $\Delta$. Note that due to its
block structure and the invertibility of the blocks, we can easily invert $J$.
The stable fixed point is located at

$\begin{pmatrix}\vec{x}^{*}\\\
\vec{z}^{*}\end{pmatrix}=J^{-1}\vec{h}=-\frac{2}{\varepsilon}\left(\mathds{1}+\frac{\Delta^{2}}{\omega_{c}^{2}}\right)^{-1}\otimes\begin{pmatrix}\mathds{1}&\frac{\Delta}{\omega_{c}}\\\
-\frac{\Delta}{\omega_{c}}&\mathds{1}\end{pmatrix}\vec{h},$ (S57)

where in the last expression, the first (diagonal) matrix is to be multiplied
with each of the four blocks (also diagonal) in the second matrix, as
indicated by the symbol $\otimes$. The fixed point solution for the components
of $\vec{\xi}$ lies at

$\xi_{\ell}=x_{\ell}^{*}\cos\omega_{c}t-z_{\ell}^{*}\sin\omega_{c}t=-\frac{\kappa_{zk}\omega_{c}}{\omega_{c}^{2}+\Delta_{\ell}^{2}}v_{k}^{[\ell]}\cos\omega_{c}t-\frac{\kappa_{zk}\Delta_{\ell}}{\omega_{c}^{2}+\Delta_{\ell}^{2}}v_{k}^{[\ell]}\sin\omega_{c}t.$
(S58)

With Eqs. (S34), (S42) and (S58), and the zeroth order approximation of
$\vartheta_{z}$, we can now evaluate Eq. (S17),

$0=Z(\omega_{s})\coloneqq
P_{z}-\alpha_{z}\omega_{s}-\frac{\kappa_{zk}^{2}}{2}\sum_{\ell=1}^{N-1}\frac{\alpha\omega_{c}\left(v_{k}^{[\ell]}\right)^{2}}{\left(\lambda^{[\ell]}-m\omega_{c}^{2}\right)^{2}+\alpha^{2}\omega_{c}^{2}}.$
(S59)

Note that we have already re-substituted $\Delta_{\ell}$ and the original
parameters by inverting Eq. (S23). Following Eq. (S20), the corotating
solitary frequency $\omega_{c}$ can be expressed by

$\omega_{c}=\left(1+\frac{\alpha_{z}}{\alpha}\frac{1}{N-1}\right)\omega_{s}.$
(S60)

We assume that $\alpha_{z}$ and $\alpha$ are of similar magnitude, hence for
very large $N$, one could approximate $\omega_{c}\approx\omega_{s}$, however,
we keep the more general form. Since the Laplacian given in Eq. (S38) and its
eigenvectors and eigenvalues depend on $\langle p_{zk}\rangle$, we need to
approximate it, in order to evaluate Eq. (S59). We assume that the dependence
is weak, and choose to use the steady-state phase angles $\varphi_{i}^{*}$
determined by Eq. (S4) instead of $\vartheta_{i}^{*}$ (for $i\in S$).
Alternatively, one could omit $\langle p_{zk}\rangle$ in the definition of
$\hat{P}_{i}$, and $\omega_{\text{sync}}$ when evaluating Eqs. (S20), (S27)
and (S38). Given an approximation of the Laplacian, Eq. (S60) closes Eq.
(S59), which can be evaluated numerically to find solutions for $\omega_{s}$.
In principle, the solutions can be used to increase the accuracy of the
approximation of the Laplacian, leading to an iterative solution algorithm.
Here, we evaluate the zeroth-order solutions of this cycle, as they are
already accurate. $L$ is then obtained by omitting column and line $z$ from
the $N\times N$ Laplacian of the total system $\mathcal{L}$, and omitting the
term in $\mathcal{L}_{kk}$ that stems from the coupling between nodes $k$ and
$z$ as well. This last steps conserves the Laplacian property. Without this
last step, due to the Poincaré separation theorem, the eigenvalues of $L$ lie
between those of $\mathcal{L}$. We suspect that something similar holds, at
least approximately, when including the correction of the diagonal entry
$\mathcal{L}_{k}k$.

The zeros of Eq. (S59) give rise to approximate solutions of the dynamics in
Eq. (S1) with $\langle\dot{\varphi}_{z}\rangle=\omega_{s}$. The linear
stability of a solution $\omega_{s}$ of $Z(\omega_{s})=0$ is given by
$\mathrm{d}Z/\mathrm{d}\omega_{s}<0$, which is equivalent to

$\displaystyle\alpha_{z}$ $\displaystyle<$
$\displaystyle-\frac{\mathrm{d}\langle
p_{zk}\rangle}{\mathrm{d}\omega_{s}}=-\frac{\mathrm{d}\langle
p_{zk}\rangle}{\mathrm{d}\omega_{c}}\times\frac{\mathrm{d}\omega_{c}}{\mathrm{d}\omega_{s}}$
(S61) $\displaystyle=$
$\displaystyle\sum_{\ell=1}^{N-1}\kappa_{zk}^{2}\alpha\left(v_{k}^{[\ell]}\right)^{2}\frac{\left(\lambda^{[\ell]}-m\omega_{c}^{2}\right)^{2}+\alpha^{2}\omega_{c}^{2}+2\omega_{c}^{2}\left[2m\left(\lambda^{[\ell]}-m\omega_{c}^{2}\right)-\alpha^{2}\omega_{c}\right]}{2\left[\left(\lambda^{[\ell]}-m\omega_{c}^{2}\right)^{2}+\alpha^{2}\omega_{c}^{2}\right]^{2}}$
$\displaystyle\quad\times\left(1+\frac{\alpha_{z}}{\alpha}\frac{1}{N-1}\right).$

Graphically, this means that the curve $\langle p_{zk}\rangle(\omega_{s})$
intersects the straight line $P_{z}-\alpha_{z}\omega_{s}$ with a slope greater
than $-\alpha_{z}$, i.e., coming from below in positive
$\omega_{s}$-direction.

###  I.3 Linear response

Here, we provide a second derivation of Eq. (S59) with a heuristic, but
intuitive method that relies on linear response, cf. reference [26]. If we
split the linearized energy flow in Eq. (S34) in its mean and oscillating
part, cf. Eq. (S16), we can approximate

$p_{zk,\text{osc}}(t)\approx\kappa_{zk}\sin\omega_{c}t,$ (S62)

in Eq. S32) for small $\vartheta_{k}$. This is valid, if the neglected
nonlinear term $\vartheta_{k}\sin\vartheta_{z}$ is only a small contribution
to the overall coupling of $\vartheta_{k}$. As noted before, the term should
not be neglected in Eq. S33), where it is the only coupling term and crucial
for nonzero $\langle p_{zk}\rangle$ in Eq. (S17).

Now with only harmonic forcing, we can apply linear response. Defining
$\chi_{i}\coloneqq\omega_{i}-\omega_{\text{sync}}=\dot{\vartheta}_{i}$, Eq.
(S3) for the synchronized cluster transforms into

$m\dot{\chi}_{i}\approx\hat{P}_{i}-\alpha\chi_{i}-\sum_{j=\in
S}\kappa_{ij}\sin\left(\vartheta_{i}-\vartheta_{j}\right)+\delta_{ik}\kappa_{kz}\sin\left(\omega_{c}t\right).$
(S63)

We linearize with respect to the $\vartheta_{i}$ around $\vartheta_{i}^{*}$,
and solve for them using linear response with the ansatz

$\vartheta_{i}=\vartheta_{i}^{*}+\kappa_{kz}\theta_{i}+\kappa_{kz}\vartheta_{i}^{(k)}(t).$
(S64)

Here, $\theta_{i}$ is the constant part, and $\vartheta_{i}^{(k)}(t)$ the
time-varying part of the response. The external perturbing signal is the
sinusoid with amplitude $\kappa_{kz}$ and frequency $\omega_{c}$. We write it
as a vector $\vec{F}^{(k)}(t)$ with components
$\delta_{ik}\kappa_{iz}\Im\left(e^{\imath\omega_{c}t}\right)$, where $\Im$ is
the imaginary part. The linearized system in vectorized form reads

$m\dot{\vec{\chi}}=\alpha\vec{\chi}-L\kappa_{kz}\vec{\vartheta}^{(k)}+\vec{F}^{(k)}(t).$
(S65)

Zhang et al. show that the linear response is valid for amplitudes of this
magnitude, see reference [26]. Following their calculation, the constant parts
of the response are the same, $\theta_{i}=\theta$ for all $i$, and fall out of
the equation, while the time-varying part of the response is

$\kappa_{kz}\vec{\vartheta}^{(k)}(t)=\Im\left[\left(L-m\omega_{c}^{2}+\imath\alpha\omega_{c}\right)^{-1}\times\vec{F}^{(k)}(t)\right].$
(S66)

An analytical expression is obtained by projecting $\vec{F}^{(k)}(t)$ onto the
orthonormal eigenbasis of $L$, i.e., inserting
$\mathds{1}_{ij}=\sum_{\ell=1}^{N-1}v_{i}^{[\ell]}\left(v_{j}^{[\ell]}\right)^{T}$,

$\kappa_{kz}\vec{\vartheta}^{(k)}(t)=\Im\left(\sum_{\ell=1}^{N-1}\frac{\vec{v}^{[\ell]}v_{k}^{[\ell]}}{\lambda^{[\ell]}-m\omega_{c}^{2}+\imath\alpha\omega_{c}}\kappa_{kz}e^{\imath\omega_{c}t}\right).$
(S67)

If we define the modal response amplitudes at node $i$ as

$a_{i}^{(k)[\ell]}\coloneqq\frac{\kappa_{kz}v_{i}^{[\ell]}v_{k}^{[\ell]}}{\sqrt{\left(\lambda^{[\ell]}-m\omega_{c}^{2}\right)^{2}+\alpha^{2}\omega_{c}^{2}}},$
(S68)

and the modal phase lags $\delta^{[\ell]}$ as

$\displaystyle\tan\delta^{[\ell]}$
$\displaystyle\coloneqq-\frac{\alpha\omega_{c}}{\lambda^{[\ell]}-m\omega_{c}^{2}},$
$\displaystyle\sin\delta^{[\ell]}$
$\displaystyle=-\frac{\alpha\omega_{c}}{\sqrt{\left(\lambda^{[\ell]}-m\omega_{c}^{2}\right)^{2}+\alpha^{2}\omega_{c}^{2}}},$
$\displaystyle\cos\delta^{[\ell]}$
$\displaystyle=\frac{\lambda^{[\ell]}-m\omega_{c}^{2}}{\sqrt{\left(\lambda^{[\ell]}-m\omega_{c}^{2}\right)^{2}+\alpha^{2}\omega_{c}^{2}}},$
(S69)

we can write the response as

$\kappa_{kz}\vartheta_{i}^{(k)}(t)=\Im\left[a_{i}^{(k)[\ell]}e^{\imath(\omega_{c}t+\delta^{[\ell]})}\right].$
(S70)

This provides us with the following interpretation. The solitary rotation with
frequency $\omega_{c}$ excites network modes given by
$\omega^{[\ell]}\approx\sqrt{\lambda^{[\ell]}}$. The nodal amplitudes depend
on the mismatch between the driving frequency $\omega_{s}$ and the
eigenfrequencies of the network, as well as the localization of the
corresponding eigenvectors and the damping.

To close the self-consistent equation, we need to determine how the excitation
of the network acts back on $\omega_{c}$. Averaging over $\dot{\omega}_{z}$,
see Eq. (S17), using Eqs. (S34) and (S70) and
$\vartheta_{z}\approx\omega_{c}t$, yields the same self-consistent equation as
the averaging approach, Eq. (S59). In terms of the response amplitudes and
phase lags, the self-consistency function takes the form

$Z(\omega_{s})=P_{z}-\alpha_{z}\omega_{s}+\frac{\kappa_{kz}}{2}\sum_{\ell=1}^{N}a_{k}^{(k)[\ell]}\sin\delta^{[\ell]}.$
(S71)

In fact, the terms in the denominators in the sum over all modes in Eq. (S59)
resemble the resonance curve of linear oscillators. It is known as
Cauchy–Lorentz distribution, Lorentz(ian) function, or Breit–Wigner
distribution.

###  I.4 Validity

The derivation is based on system parameters and topologies that result in
small $\varepsilon$ and localized eigenmodes, which is always given in our
examples. Precise but concise general statements of the validity range of our
theory would be difficult to obtain. To test the fulfillment of assumptions,
one could compare the approximation of the trajectories from Eqs. (S34), (S42)
and (S58) or statistics obtained from those with simulations. The same goes
for Eq. (S70). Note that an approximation of $\vartheta_{z}$ as in S30 can be
obtained by approximation and integration of S37. See reference [48] for some
examples. We observe that the self-consistent equation tends to have no
solutions in some regimes where it is not valid, see reference [48].

###  I.5 Solitary nodes with several neighbors

In the previous derivations, we assumed that the solitary node is of degree
one with the neighbor $k$. If the solitary node has a higher degree, we can
generalize the self-consistent equation as follows. Due to linearity, we can
simply superpose the perturbations entering the network through the links from
the solitary node to its neighbors $k\in\mathcal{N}_{z}$. This transforms Eq.
(S15) into

$m_{z}\ddot{\varphi}_{z}=P_{z}-\alpha_{z}\dot{\varphi}_{z}-\sum_{k\in\mathcal{N}_{z}}p_{zk},$
(S72)

i.e., replacing $p_{zk}$ by $p_{z}\coloneqq\sum_{k}p_{zk}$ from there on. For
the averaging approach, this corresponds to an additional sum over
$k\in\mathcal{N}_{z}$ in Eq. (S55) and thereafter. For the linear response, we
sum over several perturbation vectors in Eqs. (S63) to (S65). The result is

$0=Z(\omega_{s})\coloneqq
P_{z}-\alpha_{z}\omega_{s}-\sum_{k\in\mathcal{N}_{z}}\frac{\kappa_{zk}^{2}}{2}\sum_{\ell=1}^{N-1}\frac{\alpha\omega_{c}\left(v_{k}^{[\ell]}\right)^{2}}{\left(\lambda^{[\ell]}-m\omega_{c}^{2}\right)^{2}+\alpha^{2}\omega_{c}^{2}}.$
(S73)

Eq. (S61) also gets an additional sum over $k\in\mathcal{N}_{z}$.

###  I.6 Heterogeneous system parameters

So far, we assumed that the synchronized cluster has homogeneous inertia
$m_{i}=m$ and damping $\alpha_{i}=\alpha$ for $i\in S$, while allowing for
different parameters $m_{z}$ and $\alpha_{z}$ at the solitary node, as well as
heterogeneity in the power constants $P_{i}$ and coupling matrix $\kappa_{ij}$
for the whole network. If one allows generally heterogeneous inertia and
damping in the whole network, the self-consistent equation can still be
evaluated numerically 444Note that numerical computation of the eigenvectors
and eigenvalues of the Laplacian is necessary for complex networks, anyway.
Numerical evaluation of the zeros of the self-consistent equation is necessary
for networks as complex as the star graph, already., as outlined below.

To see this, consider Eq. (S66) and replace $m$ by
$M\coloneqq\text{diag}(m_{i})$ and $\alpha$ by
$A\coloneqq\text{diag}(\alpha_{i})$ for $i\in S$. If $M$, $A$ and $L$ do not
commute 555As Jörg Liesen pointed out in personal correspondence, for a
general connected graph Laplacian, the only commuting diagonal matrix of the
Laplacian is the identity up to a scalar multiplicative constant., they can
not be diagonalized simultaneously and the projection onto $L$’s eigenbasis of
network modes can not be applied to invert the matrix

$R^{-1}(\omega_{c})\coloneqq L-\omega_{c}^{2}M+\imath\omega_{c}A.$ (S74)

However, this can still be done numerically. Multiplication of the response
matrix $R(\omega_{c})$ with $\vec{F}^{(k)}(t)$, taking the imaginary part,
projecting out component $k$ to get $\vartheta_{k}^{(k)}(t)$, and time-
averaging Eq.(S34), which is a projection on the $\kappa_{zk}\cos\omega_{c}t$
function, all amount to taking the matrix element
$-\frac{\kappa_{kz}^{2}}{2}\Im\left(R_{kk}\right)$ for each solitary neighbor
$k\in\mathcal{N}_{z}$. Hence, the self-consistent equation can be written as

$0=Z(\omega_{s})\coloneqq
P_{z}-\alpha_{z}\omega_{s}+\sum_{k\in\mathcal{N}_{z}}\frac{\kappa_{zk}^{2}}{2}\Im\left[R_{kk}(\omega_{c})\right].$
(S75)

Note that, here, $\omega_{\text{sync}}$ does not obey Eq. (S20), but the
following generalized version. Averaging over Eq. (S18) with the ansatz
$\langle\dot{\varphi}_{i}\rangle=\omega_{\text{sync}}$ for $i\in S$ yields

$\omega_{\text{sync}}=\frac{-P_{z}+\langle p_{z}\rangle}{\sum_{i\in
S}\alpha_{i}}=\frac{-\omega_{s}\alpha_{z}}{\sum_{i\in
S}\alpha_{i}}=-\frac{\alpha_{z}}{\alpha^{*}}\frac{\omega_{s}}{N-1},$ (S76)

where $\alpha^{*}$ is the average $\alpha_{i}$ in the synchronized cluster,
and

$\omega_{c}\coloneqq\omega_{s}-\omega_{\text{sync}}=\left(1+\frac{\alpha_{z}}{\sum_{i\in
S}\alpha_{i}}\right)\omega_{s}=\left(1+\frac{\alpha_{z}}{\alpha^{*}}\frac{1}{N-1}\right)\omega_{s}.$
(S77)

We expect that the mechanism behind the solitary frequency tuning in
heterogeneous networks is still resonance with network modes, as was the case
for homogeneous networks. The modes are still the imaginary parts of the
Jacobian’s eigenvalues in Eq. (S6), however, a simple relationship to the
eigenvalues of the Laplacian like Eq. (S7) does not hold anymore. We reserve
it for future work to test Eq. (S75) for heterogeneous $M$ and $A$.

###  I.7 Several solitary nodes

The self-consistency framework presented in this work can be generalized from
1-solitary states to $l$-solitary states with $l$ nodes that oscillate around
desynchronized frequencies. If all $l$ nodes have different frequencies, they
do not interact with each other directly, because their coupling cancels out.
However, they all excite network modes, and we can superpose the respective
responses of the synchronized cluster, which has size $N-l$ now. Analogous to
the considerations in this work, we can get a set of $l$ self-consistent
equations that are only weakly coupled through $\omega_{\text{sync}}$, because
averaging over the equations for $\ddot{\varphi}_{l}$ averages out the
responses to frequencies other than the respective $\omega_{c;l}$.

If solitary nodes have identical frequencies, we enter the regime of cluster
synchronization. When averaging over the equations for $\ddot{\varphi}_{l}$,
the response of the network to all identical frequencies add up in every
respective equation. One has to take care of additional constant phase shifts
between the members of the cluster 666Note that the system is of $l$ equations
determines $\omega_{s}$ and $l-1$ relative phases, as we can fix one.. In this
work, only concerned with 1-solitaries, this could be ignored due to the
averaging process.

We expect that there are limits to the amount of perturbations and responses
that can be superposed, because the assumption that linearization works might
break down at some point. It is a task for future work to explore the
capabilities and limits of extensions of the theory of solitary states
presented here.

##  II Evaluation of the self-consistent equation

Starting from the self-consistent equation, Eq. (S59), one can apply model
reduction by simply replacing the Laplacian eigenvalues and eigenvectors with
the ones from effective models. In the following, we evaluate several
effective models and a complex network without any reduction as examples.

###  II.1 The minimal effective model

The minimal effective model is the simplest model that shows tristability
between the synchronous state and solitary states at several frequencies. The
graph of the minimal effective model can be seen in the inset in Fig. 3
(right) in the main article. Starting from a network with dense a sprout node,
$z$, and root node, $k$, we rewrite the dynamics using the local order
parameter of node $k$,

$\rho_{k}e^{\imath\Psi_{k}}\coloneqq\frac{1}{\kappa_{kz}^{*}n}\sum_{j\in\mathcal{N}_{k}}\kappa_{kj}e^{\imath\varphi_{j}},$
(S78)

where

$\kappa_{kz}^{*}n\coloneqq\sum_{j\in\mathcal{N}_{k}}\kappa_{kj},$ (S79)

so that $\kappa_{kz}^{*}$ is the mean coupling strength of the links adjacent
to node $k$ except the link $(k,z)$. This is inspired by the classical
Kuramoto order parameter, see reference [8]. Now we can rewrite Eq. (S3) for
node $k$,

$m\dot{\omega}_{k}=P_{k}-\alpha\omega_{k}-\kappa_{kz}\sin(\varphi_{k}-\omega_{s}t)-\rho_{k}\kappa_{kz}^{*}n\sin(\varphi_{k}-\Psi_{k}).$
(S80)

The complex network is then reduced to a slack node 777A slack node has
infinite mass and fixed phase. It is also called infinite-grid node. It can
absorb any amount of power instantaneously and is widely used in engineering
to balance systems. and the two nodes $z$ and $k$ by setting the inertia of
the rest of the grid to infinity. This is especially accurate for large, well-
synchronized networks. The frequencies of the nodes in the synchronized
cluster are zero and their phases constant. Without loss of generality, we set
$\Psi_{k}=0$. We approximate $\rho_{k}$ by its asymptotic temporal mean value
$\langle\rho_{k}\rangle$. Between the root and the slack node, there are $n$
lines of coupling strength $K$. Due to incoherence of the nodes in
$\mathcal{N}_{k}$, we know that $\rho_{k}<1$. This leads to an effective
coupling strength $\langle\rho_{k}\rangle\kappa_{kz}^{*}n$, which is smaller
than in the idealized case with $\kappa_{kz}^{*}n$. Therefore, we expect
resonances at smaller $\omega_{s}$ in complex networks. For the analysis of
the minimal effective model, we label the sprout node $z=1$, the root node
$k=2$, and define $K_{1}\coloneqq\kappa_{12}$, $K_{2}\coloneqq\kappa_{21}$,
and set and $K_{3}n\coloneqq\kappa_{kz}^{*}n\langle\rho_{k}\rangle$.

The equations of motion of the minimal effective model are

$\displaystyle m_{1}\ddot{\varphi}_{1}$ $\displaystyle=$ $\displaystyle
P_{1}-\alpha_{1}\dot{\varphi}_{1}-K_{1}\sin(\varphi_{1}-\varphi_{2}),$ (S81)
$\displaystyle m_{2}\ddot{\varphi}_{2}$ $\displaystyle=$ $\displaystyle
P_{2}-\alpha_{2}\dot{\varphi}_{2}-K_{2}\sin(\varphi_{2}-\varphi_{1})-K_{3}n\sin\varphi_{2}.$
(S82)

We set $m_{1}=m_{2}=1$, $P_{1}=1=-P_{2}$, $\alpha_{1}=\alpha_{2}=0.1$, and
$K=6$. In a power grid context, $K_{1}=K_{2}=K$; furthermore, $n$ gives us the
freedom to set $K_{3}=K$, too. However, this system can be analyzed in a more
general form. For the slack node, we just have $\varphi_{3}=0$ at all times.
Node 1 can have different solitary frequencies, while node 2 and 3 form the
synchronized cluster 888Node 2 can also be solitary at its natural frequency,
but we focus on node 1, for which this toy model is designed.. Due to the
presence of a slack node in the synchronized cluster, we have
$\omega_{\text{sync}}=0$, hence $\omega_{c}=\omega_{s}$. Effectively,
$\varphi_{3}$ does not appear in the equations of motion, and the synchronized
cluster only contains node 2. The steady state is given by $\varphi_{2}=0$.
The Laplacian and its eigenvectors and eigenvalues are scalars, $L=K_{3}n$,
$\vec{v}=1$, $\lambda=L$. The self-consistent equation reads

$0=P_{1}-\alpha_{1}\omega_{s}\cdot\left(1+\frac{\alpha_{2}}{2\alpha_{1}}\frac{K_{1}K_{2}\alpha_{2}}{\left(K_{3}n-m_{2}\omega_{s}^{2}\right)^{2}+\alpha_{2}^{2}\omega_{s}^{2}}\right).$
(S83)

We see that $0<|\omega_{s}|<|\Omega_{1}|$. Furthermore, the self-consistency
condition is symmetric under $\omega_{s}\to\omega_{s}$ while $P_{1}\to-P_{1}$.
The conditions can be solved for $n$ explicitly, revealing two branches that
are horizontally centered around
$\omega_{s}=\text{sign}(P_{1})\sqrt{\lambda/m_{2}}$ in the
$\omega_{s}$-$n$-plane,

$n_{\pm}=\frac{m_{2}\omega_{s}^{2}}{K_{3}}\pm\sqrt{\frac{\alpha_{2}K_{1}K_{2}}{2\alpha_{1}K_{3}^{2}}\frac{\omega_{s}}{\Omega_{1}-\omega_{s}}-\frac{\alpha_{2}^{2}\omega_{s}^{2}}{K_{2}^{2}}}.$
(S84)

The $\omega_{+}=\omega_{s}(n_{+})$ branch is stable, while the $\omega_{-}$
branch has stable and unstable parts, see Fig. 3 in the main article.

###  II.2 The star graph

A more detailed model replaces the slack node in the minimal effective model
by $n$ identical nodes with zero power generation to mimic a synchronized
power-balanced grid connected that is well-connected to the root node. The
resulting star-shaped graph with $2+n$ nodes is shown in the inset in Fig. S1
(right). The equations of motion are

$\displaystyle m_{1}\dot{\omega}_{1}$ $\displaystyle=$ $\displaystyle
P_{1}-\alpha_{1}\omega_{1}-\kappa_{12}\sin\left(\varphi_{1}-\varphi_{2}\right),$
(S85) $\displaystyle m_{2}\dot{\omega}_{2}$ $\displaystyle=$ $\displaystyle
P_{2}-\alpha_{2}\omega_{2}-\kappa_{21}\sin\left(\varphi_{2}-\varphi_{1}\right)-\sum_{j=3}^{2+n}\kappa_{2j}\sin\left(\varphi_{2}-\varphi_{j}\right),$
(S86) $\displaystyle m_{i}\dot{\omega}_{i}$ $\displaystyle=$ $\displaystyle
P_{i}-\alpha_{i}\omega_{i}-\kappa_{i2}\sin\left(\varphi_{i}-\varphi_{2}\right)\qquad\text{for
}i\in\left\\{3,\ldots,2+n\right\\}.$ (S87)

For identical $\varphi_{i\geq 3}$, the latter two equations of motion can be
reduced to

$\displaystyle m_{2}\dot{\omega}_{2}$
$\displaystyle=P_{2}-\alpha_{2}\omega_{2}-\kappa_{21}\sin\left(\varphi_{2}-\varphi_{1}\right)-\kappa_{23}n\sin\left(\varphi_{2}-\varphi_{3}\right),$
(S88) $\displaystyle m_{3}\dot{\omega}_{3}$
$\displaystyle=P_{3}-\alpha_{3}\omega_{3}-\kappa_{32}\sin\left(\varphi_{3}-\varphi_{2}\right).$
(S89)

Again, we use identical $m_{i}=m=1$, $\alpha_{i}=\alpha=0.1$ and
$\kappa_{ij}=6$ (if connected). The power constants are given by
$P_{1}=-P_{2}=1$, and $P_{3}=0$. We observe that node 3 does not have solitary
states, because $P_{3}=0$. Any initial difference in the $\varphi_{i\geq 3}$
are dissipated over time, which makes any additional links between these nodes
obsolete. For example, simulations of a network with additional all-to-all
coupling except for node $1$, show the same asymptotic behavior.

The steady state of the original network, including node 1, is given by
$P_{1}=K\sin\left(\varphi_{1}^{*}-\varphi_{2}^{*}\right)$ and
$0=P_{3}=K\sin\left(\varphi_{3}^{*}-\varphi_{2}^{*}\right)$. The synchronized
cluster has the $(1+n)\times(1+n)$ steady-state Laplacian

$L=K\begin{bmatrix}n&-1&-1&\ldots&-1\\\ -1&1&0&\ldots&0\\\
-1&0&1&\ddots&\vdots\\\ \vdots&\vdots&\ddots&\ddots&0\\\
-1&0&\ldots&0&1\end{bmatrix}.$ (S90)

Its eigenvalues and eigenvectors are

$\displaystyle\left\\{\lambda^{[\ell]}\right\\}$ $\displaystyle=$
$\displaystyle K\cdot\\{0,1,\ldots,1,n+1\\},$ (S91) $\displaystyle
v_{i}^{[1]}$ $\displaystyle=$ $\displaystyle\frac{1}{\sqrt{n+1}},\quad
v_{i}^{[n+1]}=\frac{1-\delta_{i1}(n+1)}{\sqrt{n^{2}+n}}.$ (S92)

The eigenvectors to the $n-1$ eigenvalues $K$ do not contribute to the self-
consistent equation for solitary node 1, because they are strictly zero at the
neighbor, node $2$. They can be obtained by orthogonalization of

$\tilde{v}_{i}^{[\ell]}=\frac{\delta_{i2}-\delta_{i(\ell+1)}}{\sqrt{2}}\quad\text{for
{\hbox{\ell= 2,\ldots,n}}}.$ (S93)

Note that component $i$ of the eigenvectors corresponds to $\varphi_{i+1}$ in
the original graph.

Eq. (S59) explicitly evaluates to

$\displaystyle Z(\omega_{s})$
$\displaystyle=P_{1}-\alpha_{1}\omega_{s}-\frac{K^{2}\alpha\omega_{c}}{2}\left(\frac{\frac{1}{n+1}}{\left(m\omega_{c}^{2}\right)^{2}+\alpha^{2}\omega_{c}^{2}}+\frac{\frac{n}{n+1}}{\left(K(n+1)-m\omega_{c}^{2}\right)^{2}+\alpha^{2}\omega_{c}^{2}}\right),$
(S94)

where $\omega_{c}$ is given by Eq. (S60). The solutions and their linear
stability are shown for $n=6$ in Fig. S1 (left).

Figure S1: How network modes shape the landscape of solitary attractor states
for a star-shaped graph, cf. Fig. 3 in the main article. (left) Intersections
between the straight line $P_{1}-\alpha\omega_{s}$ and the mean energy flow
$p_{s}(n=6)$ generate stable and unstable solutions according to Eq. (S59). We
use $P_{1}=1$, $\alpha=0.1$. The single non-zero peak of $p_{s}$ is centered
around the network mode close to $\sqrt{K(n+1)}$. (right) Bifurcation diagram
with bifurcation parameter $n$. The inset shows the graph. Round nodes have
$P_{i}=0$. We plot the asymptotic frequencies obtained from simulations of
random initial conditions.

The bifurcation diagram is shown in the right panel. In comparison to the
minimal effective model in Fig. 3 in the main article, we see here that Eq.
(S59) for any network with a synchronized cluster of size larger than one (not
counting infinite-grid nodes) has an additional unstable solution due to the
bulk mode at $\lambda^{[1]}=0$. Other than that, the features are similar to
the minimal effective model. However, the only contributing non-zero network
mode lies close to $\lambda^{[n+1]}=K(n+1)$.

Note that the network modes that are shown are translated into the static
frame via

$\mu^{[\ell]}\coloneqq\text{sign}(P_{1})\cdot\left(1+\frac{\alpha_{z}}{\alpha}\frac{1}{N-1}\right)^{-1}\sqrt{\frac{\lambda^{[\ell]}}{m}},$
(S95)

cf. Eq. (S60). The reason is that $\omega_{c}$ is in resonance with
$\sqrt{\lambda^{[\ell]}}$, such that $\lambda^{[\ell]}-m\omega_{c}^{2}$
becomes small and
$\text{sign}(\omega_{s})=\text{sign}(\omega_{c})=\text{sign}(P_{1})$. However,
we plot $\omega_{s}$, for which the original network modes
$\pm\sqrt{\lambda^{[\ell]}/m}$ appear shifted away from the peaks.

###  II.3 Complex synthetic power grids

Fig. S2 shows the graph of a spatially embedded synthetic power grid.

Figure S2: A synthetic power grid of $N=100$ nodes and their topological
classes according to reference [23]. Generated using a random growth model,
see reference [33]. Nodes marked with a ”$+$” (”$-$”) have $P_{i}=\pm 1$. The
grid and its three dense sprout nodes, highlighted with node numbers 74, 75
and 84, are used as examples throughout the figures of the main article and
the supplement.

Fig. S3 shows properties of the ensemble of 100 synthetic power grids with
$N=100$ nodes each.

Figure S3: Topological properties of the ensemble of 100 synthetic networks.
(left) Histogram of the number of dense sprout nodes per network. (right)
Histogram of the degrees of the respective root nodes. Each root node appears
as many times as it has dense sprout nodes as neighbors.

Fig. S4 shows all asymptotic frequencies and a histogram of the number of
desynchronized nodes for grid 1.

Figure S4: Asymptotic frequencies in a synthetic power grid. We apply a fine
mesh of 1,200 evenly spaced single-node perturbations to each of the 100 nodes
in a synthetic power grid, see Fig. S2. About 22% of the perturbations lead to
desynchronized asymptotic states. (top) The wide spectrum of possible
asymptotic mean frequencies in this highly multistable system is scattered
over the node locations. We see that the natural frequencies place approximate
bounds on the asymptotic mean frequencies. The difference to Fig. 2 in the
main article is that here, all asymptotic states are shown (not only
1-solitaries). (right) Histogram of the relative abundance of solitary degrees
(numbers of desynchronized nodes, not degrees of solitary nodes) of all the
desynchronized asymptotic states. In a single occasion, a single-node
perturbation desynchronized the whole network.

Fig. S5 shows the solutions of Eq. (S59) for $z=74$ and homogeneous parameters
$P_{i}=\pm 1$ and $\alpha_{i}=0.1$ for the solitary oscillator ($P_{74}=-1$)
and all other oscillators.

Figure S5: Solutions of Eq. (S59) for a dense sprout node ($i=74)$ in a
synthetic power grid, see Fig. S2, with standard parameters. The number of
intersections is much higher, but the simulations show qualitatively similar
results.

Fig. 4 in the main article shows the same, but for $P_{74}=-1+\Delta P$ with
$\Delta P=-1.1$ and $\alpha_{74}=0.18$. The other $P_{i}$ are shifted from
their standard values $\pm 1$ by $-\Delta P/N$ to maintain power balance. This
slight adjustment makes the analytical predictions in the figure much clearer.
The simulations confirm that the dynamics, especially the locations and
stability of the solitary states, are virtually unchanged.

Fig. S6 shows the solutions of Eq. (S59) for node $75$, which is a generator
999Note that the visibilities (alpha values) of the simulation data points are
not comparable between plots due to normalization and manual increase if
necessary.. This illustrates that these intersections plots are point-
symmetric with respect to the origin if consumers and generators are
exchanged. This corresponds to the symmetry of Eq. (S59) under $P_{z}\to-
P_{z}$ while $\omega_{s}\to-\omega_{s}$, which implies
$\omega_{c}\to-\omega_{c}$. We remark that due to other modes, the solution
can be shifted significantly from the network mode that is most localized at
the root node, cf. Fig S6.

Figure S6: Solutions of Eq. (S59) for a dense sprout node ($i=75)$ in a
synthetic power grid, see Fig. S2, with standard parameters. The landscape of
network modes can shift the main solution far from the most localized network
mode.

Fig. S7 shows a phase portrait of node 74 of grid 1, the grid in Fig. S2.

Figure S7: Phase portrait of node 74 of grid 1. Color indicates the mean
frequency after the transient following perturbations to the synchronous state
with $\Delta\varphi$ and $\Delta\omega$ to $\varphi_{74}$ and $\omega_{74}$
respectively. Yellow is the synchronous state, turquoise a resonant solitary
state, purple the natural solitary state.

Fig. S8 compares the analytical prediction of the minimal effective model with
simulation data from all 462 dense sprouts in the ensemble. Our prediction
from the model reduction that the minimal model provides an upper bound are
confirmed.

Figure S8: The minimal effective model provides an upper bound for the mean
frequencies of resonant solitary states, as seen by comparison with the
simulation results from the ensemble of synthetic power grids. A mesh of
single-node perturbations is applied to all dense sprout nodes in the
ensemble. All asymptotic mean frequencies at the dense sprout nodes are then
scattered and overlaid by the analytical solution for the effective toy model.
The branches $\omega_{\pm}$ correspond to $\omega_{s}(n_{\pm})$ in Eq. (S84).
$\sqrt{Kn}$ is the network mode at the center between the branches.
$\Omega=\pm 10$ is the natural frequency of generators or consumers,
respectively, and the asymptotic bound for the solitary frequencies.