# Submodular Input Selection for Synchronization in Kuramoto Networks

Dinuka Sahabandu1, Andrew Clark2, Linda Bushnell1, and Radha Poovendran1 1D.
Sahabandu, L. Bushnell, and R. Poovendran are with the Department of
Electrical and Computer Engineering, University of Washington, Seattle, WA
98195 USA<EMAIL_ADDRESS>Clark is with the Department of
Electrical and Computer Engineering, Worcester Polytechnic Institute,
Worcester, MA 01609 USA<EMAIL_ADDRESS>
###### Abstract

Synchronization is an essential property of engineered and natural networked
dynamical systems. The Kuramoto model of nonlinear synchronization has been
widely studied in applications including entrainment of clock cells in brain
networks and power system stability. Synchronization of Kuramoto networks has
been found to be challenging in the presence of signed couplings between
oscillators and when the network includes oscillators with heterogeneous
natural frequencies. In this paper, we study the problem of minimum-set
control input selection for synchronizing signed Kuramoto networks. We first
derive sufficient conditions for synchronization in homogeneous as well as
heterogeneous Kuramoto networks using a passivity-based framework. We then
develop a submodular algorithm for selecting a minimum set of control inputs
for a given Kuramoto network. We evaluate our approach through a numerical
study on multiple classes of graphs, including undirected, directed, and cycle
graphs.

##  I Introduction

Synchronization of coupled nonlinear oscillators has been extensively studied
in engineering, physics, and biology. Applications in engineering include
power system stability [1] and clock synchronization in wireless sensor
networks [2]. In physics, synchronization has been used to study coupled
pendulum clocks [3]. It has been employed in biology to understand the
circadian rhythms in networks of bursting neurons [4] and to investigate the
synchronous firing of cardiac pacemaker cells [5].

Kuramoto dynamics [6] is a widely accepted model used to study synchronization
in networked oscillators due to its ability to capture a variety of real-world
synchronization phenomena [7]. Under the Kuramoto model each oscillator’s
dynamics consists of two components, its own natural frequency and positively-
or negatively-weighted sinusoidal couplings with the neighboring oscillators.

Synchronization in oscillator networks has been studied under two categories:
networks with homogeneous natural frequencies (e.g., synchronization of
coupled identical pendulum clocks [3]) and networks with heterogeneous natural
frequencies (e.g., study of cardiac rhythms [5]). In Kuramoto networks with
homogeneous natural frequencies, it is known that all oscillators will
converge to the same phase angle (phase synchronization) from any initial
condition when the network is undirected and contains only positive couplings
[8]. Synchronization, however, may not be realized in homogeneous Kuramoto
dynamics on weakly connected directed graphs or under signed couplings [9].
Moreover, even with all positive couplings, heterogeneous Kuramoto networks
may fail to achieve synchronization [8]. In [10], a subset of oscillators were
chosen to drive a positively coupled Kuramoto network towards synchronization.

Achieving synchronization in any Kuramoto network requires answering the
following questions. First, given any Kuramoto network, can we find a
mechanism that drives it towards synchronization? Second, what are the
conditions to drive a given Kuramoto network towards synchronization?

In this paper we develop an analytical framework for selecting a minimum set
of control inputs with constant and identical phases and natural frequencies
to drive any given Kuramoto network towards synchronization. To the best of
our knowledge, this is the first paper to do so. We make the following
contributions.

  * •

We decompose the Kuramoto model into a negative feedback interconnection
between a nonlinear subsystem describing the local node dynamics and a linear
subsystem describing the inter-node coupling.

  * •

We derive sufficient conditions for synchronization in a broad class of
networks, including homogeneous and heterogenous oscillators, undirected and
directed graphs, and networks with signed coupling.

  * •

We prove that selecting a minimum set of control inputs to satisfy the
sufficient conditions is a submodular optimization problem, and develop a
polynomial-time approximate algorithm with provable optimality bounds.

  * •

We evaluate our approach via a numerical study.

The rest of the paper is organized as follows. Section II presents the related
work. Section III provides needed preliminaries. Section IV presents the
problem of selecting control inputs. Section V presents the sufficient
conditions for synchronization. Section VI presents a submodular algorithm for
selecting a minimum set of control inputs. Section VII provides the simulation
results. Section VIII presents conclusions.

##  II Related work

Most of the work in the literature focuses on studying the
necessary/sufficient conditions to achieve phase/frequency synchronization in
homogeneous/heterogeneous Kuramoto networks with positive couplings [8, 11,
12]. In [8, 11] order parameter and Lyapunov theorem based approach was used
to derive conditions for synchronization. In [12], conditions for frequency
synchronization in a ring of unidirectionally coupled oscillators were studied
by extending Gershgorin’s theorem. A passivity-based decomposition was
proposed in [13] to study synchronization in positively coupled, homogeneous
Kuramoto networks of undirected graphs. We generalize this approach for
deriving synchronization conditions in signed heterogeneous Kuramoto networks
under various network topologies.

Recently, synchronization in signed Kuramoto networks with
homogeneous/heterogeneous natural frequencies have received attention [14, 15,
16]. In [14] synchronization of Kuramoto oscillators with all-to-all signed
couplings were studied. In [15] conditions for frequency synchronization in
signed directed networks were derived. In [16] the conditions that ensures the
synchronization in signed acyclic directed oriented heterogeneous oscillator
networks and signed oriented cyclic homogeneous oscillator networks were
analyzed. Yet, the analytical frameworks considered in these papers do not
generalize to studying synchronization conditions in signed Kuramoto networks
with general undirected and directed graphs.

The effect of a single pacemaker (leader) on synchronization of Kuramoto
networks with positive coupling under both homogeneous and heterogeneous
natural frequencies are studied in [17, 18]. In [10] a submodular optimization
framework for selecting a set of control input nodes in order to achieve
synchronization in positively coupled Kuramoto networks was analyzed. A recent
work in [19] developed a control strategy for the problem of leader-follower
frequency synchronization in positively coupled Kuramoto networks, by
exploiting the adaptive control framework. However, these models do not
address synchronization of signed Kuramoto networks. A novel passivity-based
framework and a submodular input selection algorithm presented in this paper
enables deriving sufficient conditions for synchronization and steering a
given Kuramoto network towards synchronization.

##  III Notation and Preliminaries

In this section we first introduce the notation used in this paper. Then we
provide the background on passivity and stability results from literature that
have been used to derive the theoretical results presented in this paper.

###  III-A Notation

Let $R$ be a $m\times m$ matrix. Then $[R]_{ij}$ denotes the entry in $R$
corresponding to $i^{\text{th}}$ row and $j^{\text{th}}$ column, where
$i,j\in\\{1,\ldots,m\\}$. Let $R^{T}$ denotes the transpose of matrix $R$. The
notation $R>0$ implies that $R$ is a positive definite matrix. Let
$\lambda_{\min}(R)$ denote the minimum eigenvalue of matrix $R$. The vectors
$\mathbf{0}$ and $\mathbf{1}$ represent the all zeros and all ones vectors
with appropriate dimensions, respectively. Let $\mathbb{E}(.)$ denote the
standard expectation operator. The notations $||.||_{2}$ and $||.||_{\infty}$
denote the two-norm and infinity-norm of vectors.

###  III-B Passivity and Stability

For a dynamical system $\varSigma$ with input $u(t)$ and output $y(t)$,
passivity is defined as follows.

######  Definition 1 (Passivity, [20], Corollary 2.3)

Suppose that there exists a continuously differentiable function $V\geq 0$ and
a measurable function $d$ such that $\int_{0}^{t}{d(s)\ ds}\geq 0$ for all
$t$. Then if

$\dot{V}(t)\leq y^{T}(t)u(t)-d(t)$ (1)

for all $t$ and all $u(t)$, the system $\varSigma$ is passive.

Next we introduce the notion of $\delta$-Input Strictly Passive ($\delta$-ISP)
system in the following definition.

######  Definition 2 ($\delta$-Input Strictly Passive ($\delta$-ISP) system)

Assume there exists a continuously differentiable function $V(\cdot)\geq 0$
and a measurable function $\beta(\cdot)$ such that
$\int_{0}^{t}{\beta(s)ds}\geq 0$ for all $t\geq 0$. Then a system $\varSigma$
is said to be $\delta$-Input Strictly Passive ($\delta$-ISP) if there exists a
$\delta>0$ such that

$\dot{V}(t)\leq y^{T}(t)u(t)-\delta u^{T}(t)u(t)-\beta(t)$ (2)

for all $t$ and all $u(t)$.

The following proposition provides the conditions for $\mathcal{L}_{2}$
stability of interconnected subsystems.

######  Proposition 1 ([20], Corollary 5.3)

Consider two subsystems $H_{1}$ and $H_{2}$ in a single channel negative
feedback interconnection. If $H_{1}$ is passive and $H_{2}$ is $\delta$-ISP:

  1. 1.

Closed loop system is $\mathcal{L}_{2}$ finite gain stable.

  2. 2.

There exists a $\delta>0$ such that $\mathcal{L}_{2}$ gain of the system is
bounded above by $1/\delta$ .

The next result establishes a synchronization condition for a linear system
with time-varying system matrix $A(t)$.

######  Proposition 2 ([21], Theorem 1)

Consider the system

$\dot{x}(t)=A(t)x(t),$ (3)

where $A(t)$ is a bounded and piecewise continuous function of $t$. For every
$t$, $A(t)$ is Metzler with zero row sums. For any $\Delta>0$ and any matrix
$B$, let the $\Delta$ digraph to be defined as the digraph where an edge
$(i,j)$ exists if $B_{ij}\geq 0$.

If there is an index $k\in\\{1,\ldots,n\\}$, a threshold value $\Delta>0$ and
an interval length $\hat{T}>0$ such that for all $t\in\mathbb{R}$ the
$\Delta$-digraph associated to $\int_{t}^{t+\hat{T}}A(\tau)d\tau$ has the
property that all nodes may be reached from the node $k$, then the equilibrium
set of synchronized states is uniformly exponentially stable. Solution of Eqn.
(3) converge to a common value as $t\rightarrow\infty$.

##  IV Problem Formulation

We consider a directed Kuramoto network of $n$ oscillators indexed
$\\{1,\ldots,n\\}$ with couplings among the oscillators defined by the edge
set $E:=\\{1,\ldots,m\\}$. An ordered pair $(i,j)\in E$ denotes an edge from
oscillator $i$ to oscillator $j$ and implies that the oscillator $j$ is
influenced by the oscillator $i$. Let the underlying graph structure of the
Kurmaoto network is denoted by $\mathcal{G}=(V,E)$, where the set $V$ denotes
the set of nodes111In this paper we use the words nodes and oscillators
interchangeably. that represent the Kuramoto oscillators. Define a vector of
length $n$ by $\omega=\begin{bmatrix}\omega_{i}\end{bmatrix}^{n}_{i=1}$ where
each $\omega_{i}\in\mathbb{R}$ represents the natural frequency associated
with the $i^{\text{th}}$ oscillator. Let $N_{\text{in}}(i):=\\{j:(j,i)\in
E\\}$ denote the set of oscillators that have an incoming edge to oscillator
$i$.

Kuramoto dynamics of the $i^{\text{th}}$ oscillator is given by

$\dot{\theta}_{i}(t)=\omega_{i}-\sum_{j\in
N_{\text{in}}(i)}{{K_{ji}}\sin{(\theta_{i}(t)-\theta_{j}(t))}}$ (4)

where the coupling coefficients (i.e., edge weights) $K_{ji}$ are nonzero real
numbers that characterize the influence of oscillator $j$ on oscillator $i$’s
dynamics. We define the $n\times m$ incidence matrix $D$ by

$D_{ie}=\left\\{\begin{array}[]{ll}1,&e=(j,i)\ \mbox{for some }j\\\
-1,&e=(i,j)\ \mbox{for some }j\\\ 0,&\mbox{else}\end{array}\right.$

We define a matrix $\hat{D}$ by

$\hat{D}_{ie}=\left\\{\begin{array}[]{ll}1,&e=(j,i)\ \mbox{for some }j\\\
0,&\mbox{else}\end{array}\right.$

Next, we define $m\times m$ diagonal matrix $K$ by $K_{ee}=K_{ij}$ where
$e=(i,j)$. Under these definitions, the dynamics of the network can be written
in vector form as

$\dot{\theta}(t)=\omega-\hat{D}K\sin{(D^{T}\theta(t))}$ (5)

where $\sin{\theta}=(\sin{\theta_{i}}:i=1,\ldots,n)$.

We assume that the network operates with a set of inputs,
$S\subseteq\\{1,\ldots,n\\}$. Initial phase angles and natural frequencies
associated with input nodes are set to zero (i.e, for all $i\in S$,
$\theta_{i}(0)=0$ and $\omega_{i}=0$). Therefore, input nodes maintain a
constant phase of zero regardless of the inputs of their neighbors. In such
networks, the non-input node dynamics can be written as

$\dot{\theta}_{i}(t)=\omega_{i}-\sum_{j\in N(i)\cap
S}{{K_{ji}}\sin{\theta_{i}(t)}}-\sum_{j\in N(i)\setminus
S}{{K_{ji}}\sin{(\theta_{i}(t)-\theta_{j}(t))}},$ (6)

while $\theta_{i}(t)\equiv 0$ for all $i\in S$.

Let $E(S)$ denote the set of edges that are incoming to input nodes. In
networks with inputs, we define the $(n-|S|)\times(m-|E(S)|)$ matrix $D(S)$ as

$D(S)_{ie}=\left\\{\begin{array}[]{ll}1,&i\notin S,e=(j,i)\ \mbox{ for some
}j\\\ -1,&i\notin S,e=(i,j)\ \mbox{ for some }j\notin S\\\
0,&\mbox{else}\end{array}\right.$

The effect of $S$ on $D$ is to remove all rows of $D$ related to $S$ and all
columns of $D$ representing incoming edges to $S$. The
$(n-|S|)\times(m-|E(S)|)$ matrix $\hat{D}(S)$ is defined as

$\hat{D}(S)_{ie}=\left\\{\begin{array}[]{ll}1,&i\notin S,e=(j,i)\mbox{ for
some }j\\\ 0,&\mbox{else}\end{array}\right.$

Let $(m-|E(S)|)\times(m-|E(S)|)$ diagonal matrix $K(S)$ defined by
$K_{ee}=K_{ji}$ where $e=(j,i)$ and $i\neq S$. Furthermore, with abuse of
notation we let $\omega(S)=\begin{bmatrix}\omega_{i}\end{bmatrix}_{i\notin
S}$. Then Kuramoto model can be written as

$\dot{\theta}(t)=\omega(S)-\hat{D}(S)K(S)\sin{(D(S)^{T}\theta(t))}.$ (7)

Setting $z(t)=D(S)^{T}\theta(t)$, we arrive at the following equivalent model
that gives the dynamics in terms of the inter oscillator phase angle
differences, $z(t)$.

$\dot{z}(t)=D(S)^{T}\omega(S)-D(S)^{T}\hat{D}(S)K(S)\sin{z(t)}.$ (8)

Next we define the notion of phase and frequency synchronization in Kuramoto
networks.

######  Definition 3 (Phase synchronization)

Kuramoto network is said to achieve phase synchronization if and only if
$\lim_{t\rightarrow\infty}{\theta(t)}\rightarrow c_{p}\mathbf{1}$, where
$c_{p}\in\mathbb{R}$.

######  Definition 4 (Frequency synchronization)

Kuramoto network is said to achieve frequency synchronization if and only if
$\lim_{t\rightarrow\infty}{\dot{\theta}(t)}=c_{f}\mathbf{1}$, where
$c_{f}\in\mathbb{R}$ or equivalently,
$\lim_{t\rightarrow\infty}{\theta(t)}=\bar{c}_{f}$, where
$\bar{c}_{f}\in\mathbb{R}^{m}$.

Moreover, since the input nodes have fixed states of $0$, $c_{p}=c_{f}=0$. In
what follows, we use the word synchronization and frequency synchronization
interchangeably.

Problems Studied: Selecting a minimum set of control inputs set to ensure (i)
frequency synchronization in Kuramoto networks with homogeneous natural
frequencies (Section V-A), (ii) phase synchronization in Kuramoto networks
with homogeneous natural frequencies (Section V-B) and (iii) frequency
synchronization in Kuramoto networks with heterogeneous natural frequencies
(Section V-C).

##  V Analytical framework for selecting inputs

Below, we present sufficient conditions for synchronization in the settings
(i)-(iii) defined above.

###  V-A Frequency synchronization: Homogeneous case

Under this case, we first derive sufficient conditions for frequency
synchronization and present a framework for selecting set of control inputs to
drive Kuramoto network with homogeneous natural frequencies (i.e.,
$\omega_{i}=\omega_{j}$ for all $i,j$) to synchronization. First note that, by
switching to a rotating frame, it can be shown that we can set
$\omega_{i}\equiv 0$ for all $i\in\\{1,\ldots,n\\}$, without loss of
generality [9].

Let $M(S)\triangleq D(S)^{T}\hat{D}(S)K(S)$ and
$R(S)\triangleq\frac{M(S)+M(S)^{T}}{2}$ (Matrices $M$ and $R$ are defined
similarly by omitting the restriction to the set $S$). In order to derive
sufficient conditions for synchronization, we represent the homogeneous
Kuramoto network by a negative feedback interconnection between two subsystems
as shown in Figure 1 (with noting that $r(t)=0$ when the natural frequencies,
$\omega_{i}\equiv 0$ for all $i$). The advantage of this approach is that
$H_{1}$ is nonlinear but does not depend on the input set $S$, while the
bottom block depends on $S$ but is linear.

Figure 1: Decomposition of (8) as a single channel negative feedback
interconnection with a constant input $D^{T}(S)\omega(S)$.

In the following we present a sufficient condition for frequency
synchronization.

###### Theorem 1

A homogeneous Kuramoto network achieves frequency synchronization if $R(S)>0$.

###### Proof:

Consider the Lyapunov storage function $V(z)=\sum_{e\in
E}{(1-\cos{z_{e}(t)})}.$ This storage function satisfies

$\dot{V}(z)=\sum_{e\in
E}{[\sin{z_{e}(t)}\dot{z}_{e}(t)]}=u_{1}(t)^{T}y_{1}(t).$ (9)

Note that $y_{1}(t)=u_{2}(t)$ and $u_{1}(t)=-y_{2}(t)$. Therefore,

$\begin{split}\dot{V}(z)=&-u_{2}(t)^{T}y_{2}(t)=-u_{2}(t)^{T}M(S)u_{2}(t)\\\
=&-u_{2}(t)^{T}\frac{M(S)+M(S)^{T}}{2}u_{2}(t)\\\
=&-u_{2}(t)^{T}{R(S)}u_{2}(t)\leq 0\end{split}$

since $R(S)>0$. Hence, by LaSalle’s Invariance Principle, the state $z(t)$
converges to the set $\\{z:\dot{V}(z)=0\\}$. Since $\dot{V}(z)$ is zero if and
only if $u_{2}(t)=0$ as $R(S)>0$, this set is equivalent to
$\\{z:\sin{z}=0\\}$, which is the discrete set of points satisfying
$[z]_{e}=k_{e}\pi$ for some set of integers $\\{k_{e}:e\in E\\}$. Since the
set of zeros of $\dot{V}(z)$ is discrete, the state trajectory $z(t)$
converges to a single fixed point, and thus
$\lim_{t\rightarrow\infty}{\dot{z}(t)}=0$, implying synchronization. ∎

In order to develop a mechanism to select control inputs, we re-index the
nodes and edges in $\mathcal{G}$ as follows. We assume that the set of non-
input nodes are indexed in $\\{1,\ldots,n-|S|\\}$, and the set of input nodes
is are indexed in $\\{n-|S|+1,\ldots,n\\}$. Let $\overline{S}$ denote the set
of non-input nodes. For two distinct sets of nodes, $S_{1}$ and $S_{2}$, let
$E(S_{1},S_{2})$ denote the set of edges from the nodes in $S_{1}$ to the
nodes in $S_{2}$. We then partition the edge set as
$E=E(\overline{S},\overline{S})\cup E(S,\overline{S})\cup
E(\overline{S},S)\cup E(S,S)$ and assume that the edges are indexed in this
order. Furthermore, for a matrix $Q$ associated with graph $\mathcal{G}$ and
set of nodes $S_{3}$, $S_{4}$ and $S_{5}$, let the matrix notation
$Q^{S_{3}}_{S_{4},S_{5}}$ denote the rows and columns of $Q$ restricted to the
node set $S_{3}$ and edges from the nodes in $S_{4}$ to the nodes in $S_{5}$.

Using the new indices, we rewrite $M$ and $M(S)$ as

$M=\left(\begin{array}[]{cccc}(D_{\overline{S},\overline{S}}^{\overline{S}})^{T}\hat{D}_{\overline{S},\overline{S}}^{\overline{S}}K_{\overline{S},\overline{S}}^{\overline{S}}&(D_{\overline{S},\overline{S}}^{\overline{S}})^{T}\hat{D}_{S,\overline{S}}^{\overline{S}}K_{S,\overline{S}}^{\overline{S}}\\\
(D_{S,\overline{S}}^{\overline{S}})^{T}\hat{D}_{\overline{S},\overline{S}}^{\overline{S}}K_{\overline{S},\overline{S}}^{\overline{S}}&(D_{S,\overline{S}}^{\overline{S}})^{T}\hat{D}_{S,\overline{S}}^{\overline{S}}K_{S,\overline{S}}^{\overline{S}}\\\
(D_{\overline{S},S}^{\overline{S}})^{T}\hat{D}_{\overline{S},\overline{S}}^{\overline{S}}K_{\overline{S},\overline{S}}^{\overline{S}}&(D_{\overline{S},S}^{\overline{S}})^{T}\hat{D}_{S,\overline{S}}^{\overline{S}}K_{S,\overline{S}}^{\overline{S}}\\\
0&0\end{array}\right.$ $\left.\begin{array}[]{cccc}0&0\\\
(D_{S,\overline{S}}^{S})^{T}\hat{D}_{\overline{S},S}^{S}K_{\overline{S},S}^{S}&(D_{S,\overline{S}}^{S})^{T}\hat{D}_{S,S}^{S}K_{S,S}^{S}\\\
(D_{\overline{S},S}^{S})^{T}\hat{D}_{\overline{S},S}^{S}K_{\overline{S},S}^{S}&(D_{\overline{S},S}^{S})^{T}\hat{D}_{S,S}^{S}K_{S,S}^{S}\\\
(D_{S,S}^{S})^{T}\hat{D}_{\overline{S},S}^{S}K_{\overline{S},S}^{S}&(D_{S,S})^{T}\hat{D}_{S,S}^{S}K_{S,S}^{S}\end{array}\right)$
(10)

$M(S)=\left(\begin{array}[]{cc}(D_{\overline{S},\overline{S}}^{\overline{S}})^{T}\hat{D}_{\overline{S},\overline{S}}^{\overline{S}}K_{\overline{S},\overline{S}}^{\overline{S}}&(D_{\overline{S},\overline{S}}^{\overline{S}})^{T}\hat{D}_{S,\overline{S}}^{\overline{S}}K_{S,\overline{S}}^{\overline{S}}\\\
(D_{S,\overline{S}}^{\overline{S}})^{T}\hat{D}_{\overline{S},\overline{S}}^{\overline{S}}K_{\overline{S},\overline{S}}^{\overline{S}}&(D_{S,\overline{S}}^{\overline{S}})^{T}\hat{D}_{S,\overline{S}}K_{S,\overline{S}}^{\overline{S}}\end{array}\right)$
(11)

Theorem 2 gives a sufficient condition for synchronization via selecting set
of control input nodes.

###### Theorem 2

Removing a subset of rows and columns from matrix $R$ such that $R(S)>0$
suffices to drive a homogeneous Kuramoto network to frequency synchronization
by selecting a set of control input nodes $S$.

###### Proof:

The matrix $M(S)$ in Eqn. (11) is equal to the $2\times 2$ top-left block
submatrix of $M$ in Eq. (10). Hence $M(S)$ can be obtained from $M$ by
removing a subset of rows and columns. As consequence, the matrix $R(S)$ can
be obtained from the matrix $R$ by removing a subset of rows and columns.
Continuing this removal process until $R(S)>0$ drives the homogeneous Kuramoto
networks towards the synchronization by Theorem 1. ∎

###  V-B Phase synchronization: Homogeneous case

The following theorem gives a necessary condition for phase synchronization in
homogeneous Kuramoto networks.

###### Theorem 3

Given $\theta_{i}(t)$ for all $i\notin S$ are initialized around the
neighborhood of origin (i.e., $0$), homogeneous Kuramoto networks achieve
phase synchronization if $R(S)>0$.

###### Proof:

Consider the following linearized version of the equivalent system given in
(8) (with $\omega\equiv 0$) at the origin.

$\dot{\tilde{z}}(t)=\tilde{A}~{}\tilde{z}(t),$ (12)

where $\tilde{A}=\frac{\partial\\{-D^{T}\hat{D}K\sin
z(t)\\}}{\partial{z(t)}}|_{z(t)=0}=-M(S).$ Then linear time-invariant system
in Eqn. (12) is asymptotically stable if and only if all the real parts of the
eigenvalues of $M(S)$ are positive. $R(S)>0$ ensures this condition. Since the
system has a unique fixed point at $\tilde{z}(t)=0$, $z(t)\rightarrow 0$ when
$t\rightarrow\infty$ and hence phase synchronization is achieved. ∎

Following the similar arguments as in Theorem 2, when the network is
initialized in the neighborhood of the origin, removing a subset of rows and
columns from matrix $R$ such that $R(S)>0$ suffices to drive a homogeneous
Kuramoto network to phase synchronization via selecting a set of control input
nodes $S$.

###  V-C Frequency synchronization: Heterogeneous case

In this section we provide the sufficient conditions for attaining frequency
synchronization in heterogeneous Kuramoto networks (i.e.,
$\omega_{i}\neq\omega_{j}$ for at least one pair of $(i,j)$ such that $i\neq
j$) via selecting set of control input nodes. In what follows we write
$D^{T}\omega$ instead of $D(S)^{T}\omega(S)$ for simplicity of notation. We
discuss how to remove the dependency of $S$ in $D^{T}(S)\omega(S)$ after the
proof of Theorem 4.

Our approach is based on analyzing the $\mathcal{L}_{2}$ stability of the
system in Figure 1 via passivity. Passivity of $H_{1}$ and $\delta$-ISP of
$H_{2}$ together with the small gain theorem in Proposition 1, gives the
following bound on the $\mathcal{L}_{2}$ gain of the system.

###### Lemma 1

Let $\delta=\lambda_{min}(R(S))$. Then for all $t$ and any $z(0)$, we have
that

$\left(||\sin{z(t)}||_{L_{2}}^{[0,t]}\right)^{2}\leq\frac{1}{\delta^{2}}||D^{T}\omega||_{2}^{2}t.$

###### Proof:

Passivity of $H_{1}$ is confirmed via Eqn. (9) by choosing the same storage
function given in the proof of Theorem 1.

Next, to show $\delta$-ISP of $H_{2}$, note that

$u_{2}(t)^{T}y_{2}(t)=u_{2}(t)^{T}R(S)u_{2}(t)\geq\lambda_{min}(R(S))u_{2}(t)^{T}u_{2}(t).$

Thus by choosing $V(z)=0$, we obtain

$\dot{V}(z)=0\leq
u_{2}(t)^{T}y_{2}(t)-\lambda_{min}(R(S))u_{2}(t)^{T}u_{2}(t),$

implying $\delta$-ISP with $\delta={\lambda_{min}(R(S))}$.

By Proposition 1, we have that for any input $r(t)$,

$\frac{||\sin{z(t)}||_{L_{2}}^{[0,T]}}{||r(t)||_{L_{2}}^{[0,T]}}\leq\frac{1}{\delta}.$

Choosing $r(t)\equiv D^{T}\omega$ then yields

$||\sin{z(t)}||_{L_{2}}^{[0,T]}\leq\frac{1}{\delta}\left(\int_{0}^{T}{||D^{T}\omega||_{2}^{2}\
dt}\right)^{1/2}.$

Squaring both sides and evaluating the integral gives the desired result. ∎

Let $z_{ji}(t)=\theta_{j}(t)-\theta_{i}(t)$ denote the inter oscillator phase
angle differences between oscillator $j$ and $i$. In Lemma 2, we provide a
preliminary sufficient condition required for the frequency synchronization in
heterogeneous Kuramoto networks222Similar results have been appeared in
literature for heterogeneous Kuramoto networks with positive couplings [18,
22].

###### Lemma 2

Heterogeneous Kuramoto networks achieve frequency synchronization if the phase
angles of the oscillators are bounded for all $t$ such that
$z_{ji}(t)\in(\pi/2,3\pi/2)$ for all edges $(j,i)$ with $K_{ji}<0$ and
$z_{ji}(t)\in(-\pi/2,\pi/2)$ for all edges $(j,i)$ with $K_{ji}>0$.

###### Proof:

By taking the time derivative of Eqn. (7), we have

$\ddot{\theta}(t)=-\hat{D}(S)K_{c}(S)D(S)^{T}\dot{\theta}(t),$ (13)

where matrix $K_{c}(S)$ is a diagonal matrix with
$[K_{c}(S)]_{ee}=K_{e}\cos(z_{e}(t))$, where $e=(j,i)$ for each
$i,j\in\\{1,\ldots,n\\}$ and $i\notin S$. Note that matrix
$\hat{D}(S)K_{c}(S)D(S)^{T}$ is the weighted Laplacian matrix of the graph
$\mathcal{G}(S)=\\{V\backslash S,E\backslash E(S)\\}$ with time-varying
weights $K_{e}\cos(z_{e}(t))$.

Furthermore the following hold if each $z_{e}(t)$ satisfy the conditions given
in theorem statement: (i) for all $e\in E$ such that $K_{e}<0$ we have
$cos(z_{e}(t))<0$ which gives $K_{e}cos(z_{e}(t))>0$. (ii) for all $e\in E$
such that $K_{e}>0$, we have $\cos(z_{e}(t))>0$ which gives
$K_{e}cos(z_{e}(t))>0$.

When all the weights, $[K_{c}(S)]_{ee}$, are positive for all $t$, the matrix
$-D(S)^{T}\hat{D}(S)K_{c}(S)$ has nonnegative off diagonal entries (Metzler)
with zero row sums. Then by Proposition 2, $\ddot{\theta}(t)$ converge to a
common value as $t\rightarrow\infty$. ∎

Note that the bounds given in Lemma 2 for $z(t)$ should satisfy for all $t$.
Therefore, next we explore the initial conditions, $z_{ji}(0)$, that will
achieve these bounds for all $t$.

In what follows we assume that the initial phase angles in Eqn. (8) satisfies
the following conditions.

###### Assumption 1

The initial state $z(0)$ satisfies $z_{ij}(0)\in(\pi/2,3\pi/2)$ for all edges
$(i,j)$ with $K_{ji}<0$ and $z_{ij}(0)\in(-\pi/2,\pi/2)$ for all edges $(i,j)$
with $K_{ji}>0$.

The following theorem establishes the sufficient condition for frequency
synchronization in Kuramoto networks with heterogeneous natural frequencies.

The rest of this section is organized as follows. In Theorem 4 we present
sufficient conditions for frequency synchronization in heterogeneous Kuramoto
networks. Then we present set of additional results in Lemma 3 and Theorem 5
that are required for the proof of Theorem 4. Next, we present the proof of
Theorem 4 and discuss frequency synchronization in heterogeneous Kuramoto
networks via selecting a minimum set of control inputs. Finally, we wrap-up
this section by presenting set of graph structures that enable conditions
given in Assumption 1.

###### Theorem 4

Suppose that the input set $S$ is chosen such that $\lambda_{\min}(R(S))$ is
strictly bounded below by $||D^{T}\omega||_{2}$. If Assumption 1 is satisfied
then heterogeneous Kuramoto networks achieve frequency synchronization.

In the following lemma we show that a nonnegative continuous function is upper
bounded for all $t$ if its integral and initial conditions are upper bounded.

###### Lemma 3

Let $f$ be a nonnegative continuous function. Suppose that, whenever $z(0)$
satisfies $f(z_{0})\leq C$ for a constant $C$,

$\int_{0}^{t}{f(z(\tau))\ d\tau}\leq tC$

for all $t$. Then for any $z(0)$ satisfying $f(z(0))\leq C$, we have
$f(z(t))\leq C$ for all $t$.

###### Proof:

Suppose the result does not hold. Then there exists $z$ and $t$ such that
$f(z(t))>C$ for some time $t$ when $z(0)=z$. Since the trajectory of $z(t)$ is
continuous, $f(z(t))$ is continuous as well. Define some notations as follows:

$\displaystyle t^{\ast}$ $\displaystyle=$
$\displaystyle\inf{\\{t:f(z(t))>C\\}}$ $\displaystyle t^{\ast\ast}$
$\displaystyle=$ $\displaystyle\inf{\\{t:f(z(t))\leq C,\ t>t^{\ast}\\}}$
$\displaystyle\overline{C}$ $\displaystyle=$
$\displaystyle\sup{\\{f(z(t)):t\in[t^{\ast},t^{\ast\ast}]\\}}$
$\displaystyle\eta$ $\displaystyle=$ $\displaystyle\overline{C}-C$
$\displaystyle t_{1}$ $\displaystyle=$
$\displaystyle\inf{\\{t:f(z(t))=C+\eta/2,t_{1}>t^{\ast}\\}}$ $\displaystyle
t_{2}$ $\displaystyle=$ $\displaystyle\inf{\\{t:f(z(t))>C+\eta/2,t>t_{1}\\}}$

For any $\epsilon>0$, let $t_{0}(\epsilon)=\sup{\\{t:f(z(t))\leq
C-\epsilon,t<t^{\ast}\\}}.$ Note that $(t^{\ast}-t_{0}(\epsilon))$ goes to
zero as $\epsilon$ goes to zero. Then

$\displaystyle\int_{t_{0}(\epsilon)}^{t^{\ast\ast}}{f(z(t))\ d\tau}$
$\displaystyle=$ $\displaystyle\int_{t_{0}(\epsilon)}^{t^{\ast}}{f(z(t))\
d\tau}+\int_{t^{\ast}}^{t_{1}}{f(z(t))\ d\tau}$
$\displaystyle+\int_{t_{1}}^{t_{2}}{f(z(t))\
d\tau}+\int_{t_{2}}^{t^{\ast\ast}}{f(z(t))\ d\tau}$ $\displaystyle>$
$\displaystyle(t^{\ast}-t_{0}(\epsilon))(C-\epsilon)+C(t_{1}-t^{\ast})$
$\displaystyle+(C+\frac{\eta}{2})(t_{2}-t_{1})+C(t^{\ast\ast}-t_{2})$
$\displaystyle=$
$\displaystyle(t^{\ast\ast}-t_{0}(\epsilon))C+\frac{\eta}{2}(t_{2}-t_{1})$
$\displaystyle-\epsilon(t^{\ast}-t_{0}(\epsilon))>C(t^{\ast\ast}-t_{0}(\epsilon)).$

for $\epsilon$ sufficiently small. Let $z_{0}=z(t_{0}(\epsilon))$, so that
$f(z(t))\leq C$. Set $z(0)=z_{0}$ and choose $t=t^{\ast\ast}-t_{0}(\epsilon)$.
By time invariance, $\int_{0}^{t}{f(z(t))\ d\tau}>Ct$ is a contradiction. ∎

Theorem 5, shows that the function $\sin{z(t)}$ has lower and upper bounds,
$-1+\epsilon$ and $1-\epsilon$ for $0<\epsilon<1$ when Assumption 1 is
satisfied and the smallest eigenvalue of $R(S)$ is strictly bounded below by
$||D^{T}\omega||_{2}$ and $z(0)$.

###### Theorem 5

Suppose that the input set $S$ is chosen such that the $\lambda_{\min}$ is
strictly bounded below by $||D^{T}\omega||_{2}$. If Assumption 1 is satisfied
then $||\sin{z(t)}||_{\infty}<1$ for all $t$.

###### Proof:

Let $f(z)=||\sin{z(t)}||_{2}^{2}$. Let $\epsilon$ satisfy
$||D^{T}\omega||_{2}^{2}/\delta^{2}=1-\epsilon$. Then from Lemma 1, for all
$t$,

$\displaystyle\int_{0}^{T}{||\sin{z(t)}||_{2}^{2}\
dt}\leq\frac{1}{\delta^{2}}||D^{T}\omega||_{2}^{2}T=(1-\epsilon)T$

Then by Lemma 3,
$||\sin{z(t)}||_{\infty}^{2}\leq||\sin{z(t)}||_{2}^{2}\leq(1-\epsilon)$. ∎

In the following we provide the proof of Theorem 4 using the results from
Lemma 2 and Theorem 5.

Proof of Theorem 4: Let $e=(j,i)$ for each $i,j\in\\{1,\ldots,n\\}$ and
$i\notin S$. Then using contradiction we first prove that $z_{e}(t)$ satisfies
the bounds given in Lemma 2 if $z_{e}(0)$ satisfies the bounds in Assumption
1.

Suppose $z_{e}(t)$ does not satisfies the bounds in Lemma 2 when $z_{e}(0)$
satisfies the bounds in Assumption 1. Suppose there exists $t$ such that
$z_{e}(t)\notin[\pi/2,3\pi/2]$ for some $e$ with $K_{e}<0$. Since $z_{e}(0)$
is in this region and $z_{e}$ is continuous, there exists $t^{\ast}<t$ such
that $z_{e}(t^{\ast})\in\\{\pi/2,3\pi/2\\}$.

This, however, implies $\sin{z_{e}(t^{\ast})}\in\\{-1,1\\}$, contradicting
Theorem 5. The contradiction in the case where there exists $t$ with
$z_{e}(t)\notin[-\pi/2,\pi/2]$ for some $e$ with $K_{e}>0$ is similar. Then
the result follows from Lemma 2. $\hbox{}\hfill{\blacksquare}$

Following the similar arguments as in Theorem 2, when a heterogeneous Kuramoto
network satisfies conditions in Theorem 4, removing a subset of rows and
columns from matrix $R$ such that $\lambda_{\min}(R(S))$ is bounded below by
$||D^{T}\omega||_{2}$ suffices to attain frequency synchronization via
selecting a set of control input nodes $S$. We note that $||D^{T}\omega||_{2}$
depends on the set $S$. In order to remove this dependency, we can replace
$||D^{T}\omega||_{2}$ with an uniform upper bound independent of $S$.
Derivation of one such upper bound is given below.

###### Theorem 6

Let $\bar{\delta}$ be defined as

$\bar{\delta}=\sqrt{\sum_{(j,i)\in
E}\max\\{(\omega_{j}-\omega_{i})^{2},\omega_{i}^{2},\omega_{j}^{2}\\}}.$ (14)

Then $\max_{S}{\\{||D(S)^{T}\omega(S)||_{2}\\}}\leq\overline{\delta}$.

###### Proof:

We have that

$\begin{split}||D(S)^{T}\omega(S)||_{2}^{2}&=\sum_{(j,i)\in
E(\overline{S},\overline{S})}(\omega_{j}-\omega_{i})^{2}+\sum_{(j,i)\in
E(S,\overline{S})}\omega_{i}^{2}\\\ &\leq\sum_{(j,i)\in
E}\max\\{(\omega_{j}-\omega_{i})^{2},\omega_{i}^{2},\omega_{j}^{2}\\}.\end{split}$

The results follows by taking the square root of the both sides of the above
inequality. ∎

The set of conditions in Assumption 1 is not feasible in general for all types
of network graphs. For example, a graph with oscillators $i$ and $j$ with
differently signed couplings between $(i,j)$ and $(j,i)$ does not satisfy the
conditions in Assumption 1. Hence, next we study set of network graph
structures that satisfy the conditions given in Assumption 1.

Let $E_{p}^{l}(i,j)$ and $E_{n}^{l}(i,j)$ denote the number of positive and
negative edges in a path $l$ between the nodes $i$ and $j$ in $\mathcal{G}$
such that path length is strictly larger than one, respectively. Then the
following results in Theorem 7, Corollary 1 and Corollary 2 present set of
graphs that satisfy the conditions given in Assumption 1. Proofs of Theorem 7,
Corollary 1 and Corollary 2 are given in Appendix.

###### Theorem 7

Any undirected or directed oriented graph that fulfills one of the following
conditions for each $(i,j)\in E$ satisfies Assumption 1.

  1. 1.

If $K_{ij}>0$, then for a path $l$ such that
$(E_{p}^{l}(i,j)-E_{n}^{l}(i,j))\geq 0$ requires
$(E_{p}^{l}(i,j)-E_{n}^{l}(i,j))\mod 4\in\\{0,1\\}$ and for a path $l$ such
that $(E_{n}^{l}(i,j)-E_{p}^{l}(i,j))\geq 0$ requires
$(E_{n}^{l}(i,j)-E_{p}^{l}(i,j))\mod 4\in\\{1,3\\}$

  2. 2.

If $K_{ij}<0$, then for a path $l$ such that
$(E_{p}^{l}(i,j)-E_{n}^{l}(i,j))\geq 0$ requires
$(E_{p}^{l}(i,j)-E_{n}^{l}(i,j))\mod 4\in\\{2,3\\}$ and for a path $l$ such
that $(E_{n}^{l}(i,j)-E_{p}^{l}(i,j))\geq 0$ requires
$(E_{n}^{l}(i,j)-E_{p}^{l}(i,j))\mod 4\in\\{1,2\\}$

###### Corollary 1

Let $E_{p}$ and $E_{n}$ denote the number of positive and negative edges in a
cycle graph. Then the underlying cycle graph satisfies Assumption 1 if one the
following conditions are met.

  1. 1.

$(E_{p}-E_{n})\geq 1$ and $(E_{p}-E_{n})\mod 4\in\\{1,2\\}$

  2. 2.

$(E_{n}-E_{p})\geq 2$ and $(E_{n}-E_{p})\mod 4\in\\{2,3\\}$

###### Corollary 2

Any Tree graph will satisfy Assumption 1.

##  VI Submodular Algorithm for selecting a minimum set of control inputs

In this section we provide an algorithm for selecting a minimum set of control
inputs for phase/frequency synchronization in homogeneous Kuramoto networks
and frequency synchronization in heterogeneous Kuramoto networks.

Using the results derived in Section V, we formulate the minimum-set control
input selection problem as

$\begin{split}&\min|S|\\\
&\text{s.t.}~{}\lambda_{\min}(R(S))=\delta\end{split}$ (15)

Note that when $\delta=0$ phase and frequency synchronization is achieved in
homogeneous Kuramoto networks. On the contrary, when $\delta=\bar{\delta}$ (in
Eqn. (14)) frequency synchronization is achieved in heterogeneous Kuramoto
networks.

In what follows, we show that Problem 15 can be written as a submodular
optimization problem. The formal definition of submodularity of a set function
$f$ is given below.

###### Definition 5

A set function $f:2^{V}\rightarrow\mathbb{R}$ is submodular if and only if,
for any $S\subseteq T\subseteq V$ and any $v\notin T$,

$f(S\cup\\{v\\})-f(S)\geq f(T\cup\\{v\\})-f(T),$

where $2^{V}$ denote the set of all subsets of $V$.

The following proposition relates the problem of bounding the minimum
eigenvalue of a symmetric matrix $R$ below by a constant $\delta>0$ after
removing set of corresponding rows and columns in matrix $R$ to a function
$Q(E(S))$ which is increasing and submodular in set $E(S)$.

######  Proposition 3 ([23], Lemma 2)

Let $E:=\\{1,\ldots,m\\}$ denote the set of indices related to corresponding
rows and columns in a $m\times m$ symmetric matrix $R$. Then for any subset,
${E}(S)\subseteq E$ and define a diagonal matrix of size $m\times m$ by
$\text{{diag}}(E(S))$ with $[\text{{diag}}(E(S))]_{ii}=1$ for all $i\in E(S)$
and $[\text{{diag}}(E(S))]_{ii}=0$ for all $i\notin E(S)$. Then the following
statements are equivalent:

  1. 1.

$\lambda_{\min}(R(E\backslash{E}(S)))>\delta$.

  2. 2.

There exists a constant $\alpha>0$ such that

$\lambda_{\min}(R+\alpha\text{{diag}}(E(S)))>\delta.$

  3. 3.

If $\bm{w}$ is an $m$-dimensional Gaussian random vector with mean $0$ and
covariance matrix $I$, then

$Q(E(S)):=\mathbb{E}\big{(}\min\\{w^{T}Rw+\alpha\sum_{i\in
E(S)}w_{i}^{2},\delta\\}\big{)}=\delta.$ (16)

Using the results presented in Proposition 3, we rewrite the problem of
minimum-set control input selection as follows.

$\begin{split}&\min|S|\\\ &~{}\text{s.t.}~{}Q(E(S))=\delta\end{split}$ (17)

Let ${E}(i):=\\{e~{}|~{}e=(j,i)~{}\text{and}~{}j\in N_{\text{in}}(i)\\}$
denote the set of edges incoming to oscillator $i\in\\{1,\ldots,n\\}$ in the
Kuramoto network. Problem (15) leads to the following submodular minimum-set
control input selection algorithm. Proposition 4 provides the optimality
bounds of Algorithm VI.1.

Algorithm VI.1 Algorithm for selecting a minimum set of control inputs to
induce frequency/phase synchronization.

Input:  Symmetric Matrix $R$,  $\delta=\begin{cases}0,~{}\text{if network is
homogeneous }\\\
\bar{\delta}~{}\text{in~{}Eqn}.~{}\eqref{eqn:upper_bound},~{}\text{if network
is heterogeneous}\end{cases}$

Output:  Input node set $S$

1:$S\leftarrow\emptyset$

2:while $Q(E(S))<\delta$ do

3: $v^{*}\leftarrow\text{arg max}\\{Q(E(S)\cup E(i)):i\notin S\\}$

4: $S\leftarrow\\{v^{*}\\}$

5: Calculate $Q(E(S))$ using Eqn. (16) with $R=R(S)$

6:end while

7:return $S$

###### Proposition 4

Let $T$ denote the total number of iterations Algorithm VI.1 takes to
terminate. Then define $S_{T}$ and $S^{*}$ to be the solution returned by the
Algorithm VI.1 and optimal solution of the Problem given in (17),
respectively. The Algorithm VI.1 has the following optimality bound of
$\log\frac{\delta-\lambda_{\min}(R)}{\delta-Q(E(S_{T-1}))}$, where $S_{T-1}$
denote the set of input nodes returned by algorithm in the iteration $T-1$.

###### Proof:

Notice that Algorithm VI.1 is a greedy algorithm that solves the submodular
set covering problem in (17). From [24] we obtain the following optimality
bounds for greedy algorithm that considers the submodular function $Q$.

$|S_{T}|-|S^{*}|\leq|S^{*}|\log\frac{Q(E(S_{T}))-Q(\emptyset)}{Q(E(S_{T}))-Q(E(S_{T-1}))}$
(18)

The result follows by observing $Q(E(S_{T}))=\delta$ at the convergence of
Algorithm VI.1 and $Q(\emptyset)=\lambda_{\min}(R)$ from Eqn. (16). ∎

##  VII Numerical Study

In this section we demonstrate the performance of Algorithm VI.1 in selecting
a minimum set of control inputs for synchronizing Kuramoto networks.

In what follows, we refer to Algorithm VI.1 by submodular algorithm for the
purpose of comparing the performance of Algorithm VI.1 against two other
selection algorithms: greedy and random. In each iteration greedy algorithm
selects an oscillator $i\in\\{1,\ldots,n\\}$ such that removing rows and
columns corresponding to the oscillator $i$ maximizes the $\lambda_{\min}(R)$
until $\lambda_{\min}(R(S))>\delta$. The random algorithm selects an
oscillator uniformly at random in each iteration until
$\lambda_{min}(R(S))>\delta$. We compare the submodular, random, and greedy
algorithms with the exact optimal set of inputs, which is computed via
exhaustive search.

The simulations were performed for three types of 10 node graphs, namely,
Undirected, directed oriented and directed oriented cycle graphs. Each of the
data point in the plots of Figure 2 and Figure 3 corresponds to 100 random
realizations of each graph type considered.

In Figure 2, we choose coupling weights of each $\mathcal{G}$ randomly from
the interval $[1,5]$ and then set the number of negative edges present in
$\mathcal{G}$ according to the fraction of negative edges considered. The
results suggest that submodular algorithm outperforms the random algorithm in
all three graph types compared. The number of control inputs selected by
submodular and greedy algorithm are comparable in undirected and directed
oriented graph types. However in directed cyclic graphs we observe that
submodular algorithm outperforms both greedy and random algorithms. The
average difference between the number of control inputs chosen by subbmodular
and optimal algorithms is $0.76$.

(a)

(b)

(c)

Figure 2: For homogeneous Kuramoto networks, number of control input nodes as
a function of negative edges for 10 node graphs. We compare performance of
submodular, greedy, random and optimal algorithms with each data point
averaged over 100 realizations.

In the simulations related to the heterogeneous Kuramoto networks, we use the
parameter called Weight-Frequency (WF) parameter which captures the ratio
between inter oscillator natural frequency differences $(D^{T}\omega)$ and the
in-degree ($d_{\text{in}}(i)=\sum_{(j,i)\in E(i)}K_{ji}$) of an oscillator
$i$. WF parameter is defined as

$\text{WF}=||D^{T}\omega||_{2}/\max(d_{\text{in}}(i)).$ (19)

In Figure 3 to we consider minimum-set control input selection for
synchronization in heterogeneous Kuramoto networks. In here we set fraction of
negative edges to be $0.3$ and randomly selects the natural frequencies of
oscillators from the interval $[0,2]$. The results suggests that the
performance of submodular algorithm is very similar to the results observed in
homogeneous Kuramoto networks. The average difference between the number of
control inputs chosen by submodular and optimal algorithm is $1.39$.

(a)

(b)

(c)

Figure 3: For heterogeneous Kuramoto networks, number of control input nodes
as function of WF parameter defined in Eqn. 19 for 10 node graphs. We compare
performance of submodular, greedy, random and optimal algorithms with each
data point averaged over 100 realizations.

##  VIII Conclusions

We studied the problem of minimum-set control input selection for
synchronization of Kuramoto networks. We developed a passivity-based
analytical framework to obtain sufficient conditions for synchronization. Our
framework enables control input selection in homogeneous and heterogeneous
Kuramoto networks with signed couplings under various network typologies. We
developed a submodular optimization algorithm for selecting a minimum set of
control inputs needed to drive Kuramoto networks towards synchronization. We
simulated our algorithm on undirected, directed oriented, and directed
oriented cycle graphs. Future work includes the design of time-varying input
signals associated with a minimum set of control inputs in heterogeneous
Kuramoto networks for achieving synchronization.

## References

  * [1] F. Dörfler, M. Chertkov, and F. Bullo, “Synchronization in complex oscillator networks and smart grids,” _Proceedings of the National Academy of Sciences_ , vol. 110, no. 6, pp. 2005–2010, 2013.
  * [2] L. Schenato and G. Gamba, “A distributed consensus protocol for clock synchronization in wireless sensor network,” in _46th IEEE Conference on Decision and Control_. IEEE, 2007, pp. 2289–2294.
  * [3] M. Bennett, M. F. Schatz, H. Rockwood, and K. Wiesenfeld, “Huygens’s clocks,” _Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences_ , vol. 458, no. 2019, pp. 563–579, 2002.
  * [4] E. D. Herzog, “Neurons and networks in daily rhythms,” _Nature Reviews Neuroscience_ , vol. 8, no. 10, pp. 790–802, 2007.
  * [5] V. Torre, “A theory of synchronization of heart pace-maker cells,” _Journal of theoretical biology_ , vol. 61, no. 1, pp. 55–71, 1976.
  * [6] Y. Kuramoto, “Self-entrainment of a population of coupled non-linear oscillators,” in _International symposium on mathematical problems in theoretical physics_. Springer, 1975, pp. 420–422.
  * [7] J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, “The Kuramoto model: A simple paradigm for synchronization phenomena,” _Reviews of modern physics_ , vol. 77, no. 1, p. 137, 2005.
  * [8] A. Jadbabaie, N. Motee, and M. Barahona, “On the stability of the Kuramoto model of coupled nonlinear oscillators,” in _Proceedings of the American Control Conference_ , vol. 5. IEEE, 2004, pp. 4296–4301.
  * [9] S. H. Strogatz, “From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators,” _Physica D: Nonlinear Phenomena_ , vol. 143, no. 1-4, pp. 1–20, 2000.
  * [10] A. Clark, B. Alomair, L. Bushnell, and R. Poovendran, “Toward synchronization in networks with nonlinear dynamics: A submodular optimization framework,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 10, pp. 5055–5068, 2017.
  * [11] J. C. Bronski, L. DeVille, and M. Jip Park, “Fully synchronous solutions and the synchronization phase transition for the finite-n Kuramoto model,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 22, no. 3, p. 033133, 2012.
  * [12] J. Rogge and D. Aeyels, “Stability of phase locking in a ring of unidirectionally coupled oscillators,” _Journal of Physics A: Mathematical and General_ , vol. 37, no. 46, p. 11135, 2004.
  * [13] M. Mesbahi and M. Egerstedt, _Graph Theoretic Methods in Multiagent Networks_. Princeton University Press, 2010\.
  * [14] H. Hong and S. H. Strogatz, “Kuramoto model of coupled oscillators with positive and negative coupling parameters: an example of conformist and contrarian oscillators,” _Physical Review Letters_ , vol. 106, no. 5, p. 054102, 2011.
  * [15] A. El Ati and E. Panteley, “Synchronization of phase oscillators with attractive and repulsive interconnections,” in _2013 18th International Conference on Methods & Models in Automation & Robotics (MMAR)_. IEEE, 2013, pp. 22–27.
  * [16] R. Delabays, P. Jacquod, and F. Dorfler, “The Kuramoto model on oriented and signed graphs,” _SIAM Journal on Applied Dynamical Systems_ , vol. 18, no. 1, pp. 458–480, 2019.
  * [17] X. Li and P. Rao, “Synchronizing a weighted and weakly-connected Kuramoto-oscillator digraph with a pacemaker,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 62, no. 3, pp. 899–905, 2015.
  * [18] Y. Wang and F. J. Doyle, “Exponential synchronization rate of Kuramoto oscillators in the presence of a pacemaker,” _IEEE transactions on automatic control_ , vol. 58, no. 4, pp. 989–994, 2012.
  * [19] A. Bosso, I. A. Azzollini, and S. Baldi, “Global frequency synchronization over networks of uncertain second-order Kuramoto oscillators via distributed adaptive tracking,” in _IEEE 58th Conference on Decision and Control (CDC)_. IEEE, 2019, pp. 1031–1036.
  * [20] B. Brogliato, R. Lozano, B. Maschke, and O. Egeland, _Dissipative Systems Analysis and Control_. Springer, 2007, vol. 2.
  * [21] L. Moreau, “Stability of continuous-time distributed consensus algorithms,” in _43rd IEEE conference on decision and control (CDC)_ , vol. 4. IEEE, 2004, pp. 3998–4003.
  * [22] G. S. Schmidt, A. Papachristodoulou, U. Münz, and F. Allgöwer, “Frequency synchronization and phase agreement in Kuramoto oscillator networks with delays,” _Automatica_ , vol. 48, no. 12, pp. 3008–3017, 2012\.
  * [23] A. Clark, Q. Hou, L. Bushnell, and R. Poovendran, “Maximizing the smallest eigenvalue of a symmetric matrix: A submodular optimization approach,” _Automatica_ , vol. 95, pp. 446–454, 2018.
  * [24] L. A. Wolsey, “An analysis of the greedy algorithm for the submodular set covering problem,” _Combinatorica_ , vol. 2, no. 4, pp. 385–393, 1982.

In this section, we provide proofs of Theorem 7, Corollary 1 and Corollary 2.

Proof of Theorem 7:  In this proof we let $\mathcal{G}$ to be undirected.
Notice that solution to the set of inequalities given in Assumption 1 does
exists if all the distinct inequalities involving any two nodes $i$ and $j$ in
$\mathcal{G}$ such that $(i,j)\in E$ do not contradict with each other. Also
note that the number of distinct inequalities involving any two nodes $i$ and
$j$ in $\mathcal{G}$ is equal to the number of distinct paths between node $i$
and $j$ in $\mathcal{G}$.

First consider the condition $1)$ where $K_{ij}>0$ with
$z_{ij}(0)=\theta_{i}(0)-\theta_{j}(0)$. Then from Assumption 1 we have

$-\pi/2<z_{ij}(0)<\pi/2.$ (20)

Then by adding all the inequalities defined according to Assumption 1
corresponding to the edges in the path $l$ between node $i$ and node $j$ gives

$\begin{split}-E_{(p-n)}\pi/2<z_{ij}&<E_{(p-n)}\pi/2~{}\text{if}~{}E_{(p-n)}\geq
0\\\ E_{(n-p)}\pi/2<z_{ij}&<E_{(n-p)}3\pi/2\text{if}~{}E_{(n-p)}\geq
0,\end{split}$ (21)

where $E_{(p-n)}=(E_{p}^{l}(i,j)-E_{n}^{l}(i,j))\mod 4$ and
$E_{(n-p)}=(E_{n}^{l}(i,j)-E_{p}^{l}(i,j))\mod 4$. Hence in this case we
require $E_{(p-n)}\in\\{0,1\\}$ or $E_{(n-p)}\in\\{1,3\\}$ in order for the
solutions of the inequalities in Eqns. (20) and (21) to coincide.

Next consider the condition $2)$ where $K_{ij}<0$. Then from Assumption 1 we
have

$\pi/2<z_{ij}(0)<3\pi/2.$ (22)

Then similar to the proof of condition one we can obtain the set of
inequalities in Eqn. (21) by summing up all the inequalities corresponding to
the edges in the path $l$ between node $i$ and node $j$. Therefore, we require
$E_{(p-n)}\in\\{2,3\\}$ or $E_{(n-p)}\in\\{1,2\\}$ in order for the solutions
of the inequalities in Eqns. (21) and (22) to coincide.

Similar arguments follows when the Kuramoto network’s graph structure is a
directed oriented graph. $\hbox{}\hfill{\blacksquare}$

Proof of Corollary 1: Notice that in cycle graphs for each pair of nodes
$(i,i+1)$ for $i\in{1,\ldots,n-1}$ and pair $(n,1)$ has exactly one path whose
length is greater than one. Then the proof follows form the similar arguments
as in the proof of Theorem 7 with $l=1$. $\hbox{}\hfill{\blacksquare}$

Proof of Corollary 2:  Tree graphs has exactly one path between any pair of
nodes $(i,j)$. Hence if $(i,j)\in E$ then there exists no other path between
node $i$ and node $j$ such that path length is greater than one.
$\hbox{}\hfill{\blacksquare}$


# Cluster tomography in percolation

Helen S. Ansell  Department of Physics and Astronomy, Northwestern University,
Evanston, IL 60208  Samuel J. Frank  Department of Physics and Astronomy,
Northwestern University, Evanston, IL 60208  István A. Kovács  Department of
Physics and Astronomy, Northwestern University, Evanston, IL 60208
Northwestern Institute on Complex Systems, Northwestern University, Evanston,
IL 60208  Department of Engineering Sciences and Applied Mathematics,
Northwestern University, Evanston, IL 60208

###### Abstract

In cluster tomography, we propose measuring the number of clusters $N$
intersected by a line segment of length $\ell$ across a finite sample. As
expected, the leading order of $N(\ell)$ scales as $a\ell$, where $a$ depends
on microscopic details of the system. However, at criticality, there is often
an additional nonlinearity of the form $b\ln(\ell)$, originating from the
endpoints of the line segment. By performing large scale Monte Carlo
simulations of both 2$d$ and 3$d$ percolation, we find that $b$ is universal
and depends only on the angles encountered at the endpoints of the line
segment intersecting the sample. Our findings are further supported by
analytic arguments in 2$d$, building on results in conformal field theory.
Being broadly applicable, cluster tomography can be an efficient tool to
detect phase transitions and to characterize the corresponding universality
class in classical or quantum systems with a relevant cluster structure.

Cluster formation is a prevalent feature of complex systems, including but not
restricted to magnetic domains [1], motility-induced phase separation [2, 3,
4], bacteria swarming [5], cell migration [6, 7] and the collective motion of
animal groups [8, 9, 10, 11]. As phase transitions often lead to changes in
the characteristics of the emerging clusters, phase transitions can — in
principle — be detected and studied through cluster statistics. As a
realization, let us consider a scenario where a probe is shot through a
complex system consisting of clusters and ask the following question: _How
many clusters are encountered by the probe during the measurement?_ Here we
propose that such measures of “cluster tomography” yield simple and efficient
methods for locating critical points in complex systems, both experimentally
and computationally, while also providing deep universal information on the
nature of the observed transitions.

In this letter, we first address all ten configuration types of cluster
tomography in a two-dimensional (2$d$) square system, as illustrated in Fig.
1. We focus on Bernoulli percolation on a square lattice, although the main
results are expected to apply more generally [12], as also illustrated for
3$d$ percolation. Percolation is a fundamental model of critical phenomena,
where sites or bonds of a lattice are independently occupied with probability
$p$. Collective behavior emerges as we are interested in the statistics of
clusters of connected sites [13]. In dimensions $d\geq 2$ there is a scale-
invariant critical point at $p=p_{c}$ at which the system is also believed to
be conformally invariant [14, *tsai2013]. In the most studied 2$d$ case, there
are many exact results available [16, 17, *SmirnovWerner2001], including
correlation functions [19], crossing probabilities [20] and critical
exponents.

Figure 1:  (a) Summary of all line segment types in a 2$d$ square system with
free boundaries. Endpoints of the lines are colored according to the
classification in (b), which are bulk (green), surface (orange), and
traversing (blue). (c-e) Schematics of contributing clusters and numerical
techniques used to access the critical corner contribution of different line
segments for site percolation, illustrated for line types 2, 6 and 9. Dashed
black lines represent PBC, solid lines represent FBC, and solid yellow lines
represent the measured line segment. Sites belonging to clusters that are
counted by the line segment are colored while the remaining occupied sites are
gray.

Motivated by questions on quantum entanglement in disordered systems [21, 22,
12], previous studies have considered the number of clusters (magnetic
domains), $N_{\Gamma}$ intersecting a contour, $\Gamma$. Cluster tomography
corresponds to the simplest case of skeletal entanglement [23] when $\Gamma$
is a line segment. For select line configurations of length $\ell$, it was
found that at criticality

$N(\ell)=a\ell+b\ln(\ell)+O(1)\;,$ (1)

in both 2$d$ and 3$d$ [21, 22, 24, 25]. The first term represents the expected
“area law” scaling of the number of clusters with the length of the line, with
$a$ being the non-universal linear cluster density.

The term $b\ln(\ell)$ is the “corner contribution” that emerges due to the
geometric singularities at the endpoints of the line segment, akin to a notion
of “geometric susceptibility” [26]. This term is present only at criticality
and is an elegant measure of the concavity of the underlying cluster geometry
that encapsulates universal information from correlations in the cluster shape
at all orders. For sufficiently long lines — i.e., $\ell/L=O(1)$, with $L$
being the system size — each endpoint is expected to contribute to the
“cluster count exponent” $b$ independently. For a given configuration, $b$ is
expected to be universal, that is independent of the microscopic details of
the model at criticality. The value of $b$ depends on the endpoint
configurations of the line segment, through the angle(s) appearing at each
endpoint. If the endpoint touches the system boundary, the angle $\omega$ of
the boundary at the endpoint also plays a role, see Fig. 1b. For line types
1-5 this means that the value of $b$ is universal, while for lines 6-10 $b$ is
a universal function of the characteristic angle $\gamma$.

The contribution of each endpoint depends on its topology, with three distinct
cases in 2$d$: bulk $b$, surface $s$ and traversing $t$, as illustrated in
Fig. 1b. Bulk endpoints have $\gamma=2\pi$ and occur when the endpoint is not
on the surface, while surface endpoints, for which the corner contribution
depends on $\omega$, occur when the line segment is along the surface.
Traversing endpoints occur when the endpoint is on the surface and the line
segment crosses the bulk of the system. Supported by analytic arguments from
conformal field theory (see Appendix A for details), the contributions of the
endpoints consist of linear combinations of a bulk term $B$ and a surface term
$S$, depending on the configuration. The bulk term

$B(\gamma)=\frac{5\sqrt{3}}{96\pi}\Big{(}\frac{\gamma}{\pi}-\frac{\pi}{\gamma}\Big{)}$
(2)

follows from the celebrated Cardy-Peschel formula [27], which appears whenever
the line segment passes through the bulk of the system. The surface term

$S(\gamma)=\frac{\sqrt{3}}{8\pi}\frac{\pi}{\gamma}$ (3)

appears whenever there is a change of boundary condition at the endpoint [28,
29]. Such a change of boundary condition occurs when the line segment along
which clusters are counted (fixed boundary) meets the surface (free boundary)
at an angle, which occurs for type $s$ and $t$ endpoints. Each surface
endpoint contributes $S(\omega)$ to $b$; while bulk endpoints each contribute
$B(2\pi)$. Note that usually a bulk corner contribution would be the sum of
$B(\gamma)$ and that of the conjugate angle $B(2\pi-\gamma)$ [22, 12].
However, for $\gamma=2\pi$, the contribution is just $B(2\pi)$, as the
singular $B(0)$ is not sampled by the measurement [22, 12].

Traversing endpoints are more complicated, as multiple angles contribute. The
naïve expectation would account for the bulk contributions from the angle
$\gamma$ and the complementary angle $\omega-\gamma$ between the measuring
line and the surface, as well as surface contributions from these angles due
to the boundary condition change. However, the sum of the terms
$B(\gamma)+S(\gamma)+B(\omega-\gamma)+S(\omega-\gamma)$ would correspond to a
different configuration, where clusters are also counted along the free
surface lines. Such counting along two lines at an angle $\omega$ would
correspond to $B(\omega)$ on its own, as there is no boundary condition change
(and no exterior angle). Therefore, the naïve formula is expected to be the
sum of $B(\omega)$ and the contribution of the traversing endpoint, leading to
each traversing endpoint contributing

${T}(\gamma,\omega)=B(\gamma)+S(\gamma)+B(\omega-\gamma)+S(\omega-\gamma)-B(\omega)\;.$
(4)

Note that since $B(\pi)=0$, here the $-B(\omega)$ contribution to
$T(\gamma,\omega)$ can only be detected if the endpoint is in the square
corner ($\omega=\pi/2$), corresponding to a type $t_{c}$ endpoint. Once again,
$\gamma=0$ requires special care, as $T(0,\omega)$ has unmeasured singular
contributions, so $S(\omega)$ should be used directly in this case. These
arguments lead to our predicted expressions for the cluster count exponent
$b_{ij}$ listed in Table 1, where $i,j\in\\{b,s,t\\}$ denotes the endpoint
type and an additional subscript $c$ is used on $s$ and $t$ to indicate if the
endpoint is in the corner of the square system.

These predictions for 2$d$ percolation need to be tested numerically.
Moreover, in other systems lacking detailed analytic results, numerical
measurements are the only option to perform cluster tomography. It is,
however, not immediately obvious that the universal logarithmic corner
contribution can be measured to high-precision. In practice, there are
multiple difficulties to overcome, including the fact that the linear term
$a\ell$ is much larger than the nonlinear term $b\ln(\ell)$, large statistical
noise in the overall measurement $\sim O(\sqrt{\ell})$, and potentially strong
finite-size corrections to the linear term. Strikingly, the area law term can
be canceled out _exactly_ together with its statistical error and finite-size
corrections, leading to precise measurements of the corner term even in
relatively small systems.

For line types 1 and 2, the corner contribution can be precisely measured
using a _geometric method_ [22, 24], as illustrated in Fig. 1c for a line
segment of length $\ell=L/2$ on a free surface. The difference between the
total number of clusters on two contributing segments of length $\ell=L/2$ and
the number of clusters on a full periodic line of length $L$ (without any
endpoints) gives twice the corner contribution of the line segment. As the
same sites are visited in both cases, the linear area law term cancels out and
the only remaining contribution must come from the endpoints. The corner
contribution is therefore simply half the number of shared clusters between
the two line segments. This approach has previously been used to numerically
calculate $b_{ss}$ in 2$d$, as well as $b_{bb}$ in 2$d$ and 3$d$ [22, 24],
using the same logic for a line segment in the bulk (Fig. 1a, line 1) with
periodic boundary conditions (PBC). These 2$d$ results are consistent with
predictions from conformal field theory [22], and the 2$d$ surface result has
been proven rigorously on the triangular lattice [30].

As the applicability of the geometric method alone is limited to these two
cases, here we propose a _boundary changing method_ to measure the corner
contribution of the remaining line types with high precision. This method is
expected to be applicable for cluster-based systems with short-range
interactions, and is illustrated in Fig. 1d for a line segment with traversing
endpoints on opposite sides of the system. By changing the free boundary
conditions (FBC) into PBC, the same line now forms a closed loop with no
endpoints. The area law term can be canceled by taking the difference between
the number of clusters that intersect the line segment with FBC and PBC. The
boundary changing method can be used alone, or in conjunction with the
geometric method (as illustrated in Fig. 1e), to numerically determine the
corner contributions of line types 3-10, including the angle dependence of
lines 6-10, as indicated in Table 1.

Table 1: Overview of results for 2$d$ percolation on a square lattice. Analytic and numerical values are given for the least acute angle studied $\gamma(n_{min})$ in cases where $b$ depends on $\gamma$. Parentheses by measured values give the measurement error in the last digit. “G” and “BC” respectively refer to the geometric and boundary changing methods, * indicates cancellation of lines of type 6 and/or 9 is required to determine $b$ to high precision numerically.  Line |  Endpoints111(′) for line 7 indicates endpoints are on adjacent edges | Fig. | Formula | $\gamma(n_{min})$ |  $b$ expected |  $b$ measured | Technique
---|---|---|---|---|---|---|---
1 | $bb$ |  | $2B(2\pi)$ | - |  $\frac{5\sqrt{3}}{32\pi}\approx 0.086$222Reference [22] |  0.086(1)22footnotemark: 2 | G
2 | $ss$ |  | $2S(\pi)$ | - |  $\frac{\sqrt{3}}{4\pi}\approx 0.138$333References [31, 21, 30, 22] |  0.138(3)444Reference [32] | G
3 | $s_{c}s_{c}$ |  | $2S(\frac{\pi}{2})$ | - | $\frac{\sqrt{3}}{2\pi}\approx 0.276$ | 0.27(2) | BC
4 | $ss_{c}$ |  | $S(\pi)+S(\frac{\pi}{2})$ | - | $\frac{3\sqrt{3}}{8\pi}\approx 0.207$ | 0.21(1) | G & BC
5 | $t_{c}t_{c}$ |  | $2T\left(\frac{\pi}{4},\frac{\pi}{2}\right)$ | - | $\frac{11\sqrt{3}}{8\pi}\approx 0.758$ | 0.78(2) | BC
6 | $tt$ |  | $2T(\gamma,\pi)$ | $\frac{\pi}{2}$ | $\frac{11\sqrt{3}}{16\pi}\approx 0.379$ | 0.38(1) | BC
7 | $tt^{\prime}$ |  | $T(\gamma,\pi)+T\left(\frac{\pi}{2}-\gamma,\pi\right)$ | $\frac{\pi}{4}$ | $\frac{127\sqrt{3}}{144\pi}\approx 0.486$ | 0.48(2) | BC*
8 | $tt_{c}$ |  | $T\left(\gamma,\frac{\pi}{2}\right)+T\left(\frac{\pi}{2}-\gamma,\pi\right)$ | $\arctan{\left(\frac{1}{2}\right)}$ | $\approx 0.649$ | 0.66(3) | BC*
9 | $tb$ |  | $T(\gamma,\pi)+B(2\pi)$ | $\frac{\pi}{2}$ | $\frac{27\sqrt{3}}{64\pi}\approx 0.233$ | 0.24(1) | G & BC
10 | $t_{c}b$ |  | $T\left(\gamma,\frac{\pi}{2}\right)+B(2\pi)$ | $\frac{\pi}{4}$ | $\frac{49\sqrt{3}}{64\pi}\approx 0.422$ | 0.43(1) | G & BC*

In our large-scale Monte Carlo calculations, we investigate cluster tomography
for critical site and bond percolation on the simple cubic lattice in 2$d$ and
3$d$. The respective critical occupancies for site and percolation are
$p_{c}=0.592746$ and $p_{c}=0.5$ in 2$d$ [13], and $p_{c}=0.311608$ and
$p_{c}=0.2488126$ in 3$d$ [33, *perc_cr2, 35]. We studied systems up to linear
size $L=512$ in 2$d$ and $L=256$ in 3$d$, with at least $10\,000$ samples in
each case. For each line configuration, we calculate finite-size estimates of
$b(L)$ through two-point fits at sizes $L$ and $L/2$. We then estimate $b$
using a linear extrapolation of the largest sizes against $1/L$, as shown in
Fig. 2 for a line segment spanning the entire system (line 6) in 2$d$ and
3$d$. The extrapolated $b$ values for site and bond percolation agree within
the error, indicating universality. Therefore, all quoted numerical $b$ values
are averaged for site and bond percolation.

In a 2$d$ square system, line types 1–5 each have universal $b$ values with no
angle-dependence. Numerical estimates for lines 3–5 are listed in Table 1
along with known results for lines 1–2 [31, 21, 22, 30, 32]. The angle
dependence of lines 6–10 can be explored to high precision for angles
$\gamma=\arctan(1/n)$ for integer $n_{min}\leq n\leq 20$, where $n_{min}\to 0$
for lines 6 and 9, $n_{min}=1$ for lines 7 and 10, and $n_{min}=2$ for line 8.
For line 6, we can directly access the $\gamma$ dependence of $b$ to high
precision using the boundary changing method by starting with FBC on one
opposite pair of system edges and PBC on the other. This approach can be used
in conjunction with the geometric method for line 9. For lines 7, 8 and 10,
the initial configuration must have FBC on all sides. In these cases, $b$
values can be determined numerically to high precision by incorporating lines
of type 6 and 9 into the initial configuration so that a closed loop forms
after applying PBC, then subtracting their contribution from the count to
determine the corner contribution of interest, see Fig. S1 of the Supplemental
Material [36]. Numerical results for lines 6-10 are plotted in Fig. 3a-b,
along with the predicted analytic form of the $\gamma$ dependence, and
numerical values for $\gamma(n_{min})$ for each line type are listed in Table
1. In all cases, we observe a good agreement between numerical values and
analytic predictions.

Due to the $-B(\omega)$ term in $T(\gamma,\omega)$ contributing only for type
$t_{c}$ endpoints, we further test for the presence of this term by examining
$b_{t_{c}t_{c}}$ for line 5 as the square system is sheared with acute angle
$\omega$. Analytically, we expect $b_{t_{c}t_{c}}=2T(\gamma,\omega)$ with
$\tan(\gamma)=\tan{\omega}/(1+\tan{\omega})$. Numerically, this can be tested
for $\tan{\omega}=1/n$ with $\tan{\gamma}=1/(n+1)$ in a system of size $L$ by
$L+n$. The results, plotted in Fig. 3c, show good agreement between
predictions and numerical values. As another test of our analytic predictions,
we also consider the cluster count around the entire perimeter of the sheared
square system. Since there are no bulk corners and no boundary changes in this
case, the corner contribution is expected to be zero, as demonstrated
numerically in the inset of Fig. 3c.

Figure 2: Extrapolation of the cluster count exponent $b_{tt}$ for angles
$\gamma=\arctan{(1/n)}$ for a line segment spanning the entire system (Fig. 1,
line 6). As an illustration of universality, the results are shown for both
site (  ) and bond (  ) percolation in (a) 2$d$ and (b) 3$d$.  Figure 3: (a-b)
Angle dependence of the cluster count exponent $b$ in 2$d$ for lines 6-10 in
Fig. 1a. Numerically determined points and curves showing the analytic
predictions given in Table 1 are plotted for (a) $b_{tt}$, $b_{tt^{\prime}}$
and $b_{tb}$, and (b) $b_{tt_{c}}$ and $b_{t_{c}b}$. (c) $b_{t_{c}t_{c}}$
(line 5) in a sheared square system with acute angle $\omega$. (inset) The
value of $b$ for a contour around the perimeter of a sheared square is zero.
(d) Angle dependence of $b$ in 3$d$ for $b_{tt}^{(3d)}$ and $b_{tb}^{(3d)}$.
Plotted lines show a fit of the form given in Eq. (6). (inset) Estimates of
$b$ along a partial line segment on the surface of a cube (orange up
triangles), a full line segment through the center of the face with
$\gamma=\pi/2$ (purple down triangles), and a full line along the edge of a
cube (green pentagons) all tend to zero at large $L$.

In a 3$d$ cubic system, there are 23 types of line segment configurations
depending on whether the endpoints are touching faces, edges, corners, or the
bulk, as illustrated in Fig. S2 of the Supplemental Material [36]. Of these,
only the case where both endpoints of the line segment are in the bulk has
been studied numerically, giving $b_{bb}^{(3d)}=0.130(3)$ [24]. Although there
are no known analytic predictions in 3$d$, endpoints are still expected to
contribute to $b$ independently. Therefore relations such as

$b_{tb}=(b_{tt}+b_{bb})/2$ (5)

will still hold true in 3$d$ percolation and other clustered systems. We
explore this relation by examining the angle dependence of a line segment
spanning the cube and ending on opposite faces ($b_{tt}^{(3d)}$), and the case
where the line segment starts on a face and ends in the bulk
($b_{tb}^{(3d)}$), as shown in Fig. 3d. Motivated by the expressions for
$b_{tt}$ and $b_{tb}$ in 2$d$, we fit a line of the form

$b^{(3d)}(\gamma)=c_{1}\left(\frac{\pi}{\gamma}+\frac{\pi}{\pi-\gamma}\right)+c_{2},$
(6)

which gives $c_{1}\sim 0.06$ and $c_{2}\sim 0.03$ for $b_{tt}^{(3d)}$ and
$c_{1}\sim 0.03$ and $c_{2}\sim 0.08$ for $b_{tb}^{(3d)}$. Equation (5) is
satisfied within the numerical uncertainty on the data points. However, some
line types might not to lead to a finite $b$ in 3$d$. In particular, the
critical 3d percolation occupancy is below that in 2d at which the surface
alone would lead to a logarithmic corner contribution. We investigate three
cases of a line segment on a surface of the cube: a partial line segment in a
face, a full line segment along an edge, and a full line segment through the
center of the face parallel to an edge. As shown in the inset of Fig. 3d,
these cases suggest that the cluster count exponent vanishes in 3$d$ for lines
fully on the surface (types 1-10 in Fig. S2 of the Supplemental Material
[36]). In Appendix B and Fig. S3 of the Supplemental Material [36], we discuss
how this finding is related to the shape of the clusters.

The presence of a corner contribution term at criticality provides a simple
and effective way of pinpointing the critical point and the corresponding
universality class of a clustered system. Cluster tomography with traversing
endpoints is especially promising in this regard as it leads to the strongest
signals in both 2$d$ and 3$d$, even more so at sharper angles. Although our
results are presented for percolation, which is the $Q\to 1$ limit of the
$Q$-state Potts model, similar logarithmic terms are expected to exist for
other (not necessary integer) values of $Q$, such as for the Ising model with
$Q=2$, as long as the transition is second order [12]. Techniques from
conformal field theory are expected to give the correct analytic formulas also
in these cases, at least for Fortuin-Kasteleyn clusters [37], as explored in
2$d$ for lines 1 and 2 in Ref. [12].

Our results readily provide the universal term of skeletal entanglement in the
critical bond-diluted transverse-field Ising model [38, 24], with 8 line types
studied here for the first time in 2$d$. In this application, the natural
quantity to study is the number of _crossed_ clusters, where only those
clusters contribute that are not fully contained by the line. Our numerical
findings confirm the expectation that this difference is irrelevant at the
critical point, yielding the same asymptotic behavior in the quantum model.

It would also be interesting to explore cluster tomography for percolation in
higher dimensions, as well as in the mean-field limit, although there are
clearly no corner contributions on the Bethe lattice. In general, for higher
dimensional systems, or systems lacking conformal invariance, the angle-
dependence could be nontrivial, meaning that different configurations of
cluster tomography could unveil distinct universal information.

In some cases, including experimental realizations, it may not be possible to
implement changes in the boundary conditions, meaning that one or both
numerical cancellation techniques may not apply. However, in many such cases
the nonlinear correction may still be measured to high precision while in
other cases the nonlinear term can be accessed through comparison of direct
measurements of $N(\ell)$ at different system sizes, leading to cancellation
of the linear term on average. In Fig. S4 of the Supplemental Material [36] we
demonstrate that the presence of a nonlinear corner contribution and the
corresponding $b$ exponent can be determined even without the advanced
cancellation techniques, suggesting applicability of this methodology in a
broad range of cluster-based systems.

We thank W. Witczak-Krempa for helpful discussions. We would like to
acknowledge the WCAS Summer Grant Award from the Weinberg College Baker
Program in Undergraduate Research at Northwestern University. This work was
supported by the National Science Foundation under Grant No. PHY-2310706 of
the QIS program in the Division of Physics. This research was supported in
part through the computational resources and staff contributions provided for
the Quest high performance computing facility at Northwestern University which
is jointly supported by the Office of the Provost, the Office for Research,
and Northwestern University Information Technology.

##  Appendix A Analyic arguments

Bond percolation can be considered as the $Q\to 1$ limit of the $Q$-state
Potts model [39], defined on a lattice with sites $i=1,2,\dots,n$ and $m$
nearest neighbor bonds. The partition sum of the Potts model is given by

$Z(Q)=\sum_{s}\prod_{\left\langle
ij\right\rangle}\exp\left(K\delta_{s_{i},s_{j}}\right)\;,$ (7)

where $\delta_{s_{i},s_{j}}$ is the Kronecker symbol and $K$ is the reduced
coupling, which is the ratio of the pair interaction and the temperature.
Using the identity
$\exp\left(K\delta_{s_{i},s_{j}}\right)=1+\frac{p}{1-p}\delta_{s_{i},s_{j}}$
with $p=1-e^{-K}$, the sum of products in $Z$ is written in terms of the so-
called Fortuin-Kasteleyn clusters [37], denoted by $F$. In $F$ an edge of the
lattice $i,j$ is occupied if a factor $\frac{p}{1-p}\delta_{s_{i},s_{j}}$ is
present and if the spins exist in the same state within any connected cluster.
Up to an irrelevant prefactor, this leads to

$Z(Q)\sim\sum_{F}Q^{N_{\mathrm{tot}}(F)}p^{M(F)}{(1-p)}^{m-M(F)}\;.$ (8)

For a given element of $F$ there are $N_{\mathrm{tot}}(F)\leq n$ connected
components and $M(F)\leq m$ occupied bonds. The mean total number of clusters
is then given by

$\left\langle N_{\mathrm{tot}}\right\rangle=Q\dfrac{\partial\ln Z(Q)}{\partial
Q}\;.$ (9)

If we fix all spins on a one-dimensional contour $\Gamma$ (in state $1$, say),
but leave the couplings unchanged, this relation is modified as

$\left\langle N_{\mathrm{tot}}-N_{\Gamma}\right\rangle=Q\dfrac{\partial\ln
Z_{\Gamma}(Q)}{\partial Q}\;.$ (10)

At the critical point, $p=p_{c}$, we can write [22]

$\ln Z(Q)-\ln Z_{\Gamma}(Q)\sim
f_{\mathrm{e}}(Q)L_{\Gamma}+f_{\mathrm{c}}(Q),$ (11)

where $L_{\Gamma}$ is the linear extension of $\Gamma$, $f_{\mathrm{e}}$ is
the edge free-energy density, which is a non-universal quantity, and the term
$f_{\mathrm{c}}(Q)$ originates from the corners (or endpoints) of $\Gamma$ and
is known as the _corner contribution_ to the free-energy, which is expected to
be universal. The corner contribution to the free energy can be expressed as
[27, 28, 29]

$f_{\mathrm{c}}(Q)=\bigg{(}c(Q)A_{\Gamma}+\frac{\pi}{\omega_{\Gamma}}h(Q)\bigg{)}\ln
L_{\Gamma}\;.$ (12)

The first term is the Cardy-Peschel term [27], which receives contributions
from each $\gamma_{k}$ corner of $\Gamma$, considering both the interior and
exterior sides of the contour. Here $c(Q)$ is the central charge of the
$Q$-state Potts model, and $A_{\Gamma}$ is a purely geometric factor that is
the same for all values of $Q$ and emerges from the angles $\gamma_{k}$ as

$A_{\Gamma}=\sum_{k}A_{\gamma_{k}}=\sum_{k}\frac{1}{24}\left(\frac{\gamma_{k}}{\pi}-\frac{\pi}{\gamma_{k}}\right)\;.$
(13)

The second term in Eq. (12) is present each time there are different boundary
conditions along the contour on either side of the corner, such as a change
between fixed and free boundary conditions when a measuring contour touches a
free surface. The term $h(Q)$ is the scaling dimension of the boundary
condition changing operator and $\omega_{\Gamma}$ includes the corner angle(s)
at the location of the boundary condition change.

In this letter, we focus on a line segment of length $L_{\Gamma}=\ell$ and
determine the corner contribution to the cluster number count. Using Eqs. (9)
and (10) we can write

$N(\ell)=Q\dfrac{df_{\mathrm{e}}(Q)}{dQ}\ell+Q\dfrac{df_{\mathrm{c}}(Q)}{dQ}\equiv
a\ell+b\ln{\ell}\;,$ (14)

with $b$ being the cluster count exponent. Percolation corresponds to $Q\to
1$, so using the parametrizations $c(k)=1-6/[k(k+1)]$ and
$\sqrt{Q}=2\cos(\frac{\pi}{k+1})$ [40] gives $c^{\prime}(1)=5\sqrt{3}/(4\pi)$,
where prime notation indicates a derivative with respect to $Q$. Additionally,
changing the boundary condition from fixed along the line segment to free
along a boundary gives $h^{\prime}(1)=\sqrt{3}/(8\pi)$ [31]. We conjecture
that $b$ receives contributions from linear combinations of
$B(\gamma)=c^{\prime}(1)A_{\gamma_{k}}$ and
$S(\gamma)=h^{\prime}(1)\pi/\gamma$, corresponding to Eqs. (2) and (3) of the
main text, depending on the endpoint types (surface, bulk, or traversing) of
the line segment of interest, as discussed in the main text.

##  Appendix B Relation between cluster shape and $b$

The gap-size statistics $n(s)$ considers the distance $s$ between successive
occurrences of a given cluster along a line segment through a clustered
system, quantifying the frequency with which a gap of size $s$ occurs. This
quantity provides a measure of the concavity of the cluster, and captures
pairwise and higher order correlations in the structure. For line types 1 and
2, the gap-size statistics is related to the corner contribution along a line
segment of length $\ell\leq L/2$ with PBC through [41]

$C(\ell)=\frac{1}{L}\sum_{i=1}^{\ell}\sum_{s=i}^{L/2}n(s).$ (15)

In all previously studied cases, at criticality $n(s)\sim s^{-\zeta}$ with
$\zeta=2$ [41, 32, 24]. After approximating the sum as an integral, this
results in the logarithmic corner contribution to the cluster number count in
Eq. (1).

Measures of cluster tomography on the surface of a 3$d$ cube at criticality
suggest $b=0$, indicating that $n(s)$ might decay faster. Exploring $n(s)$ of
line segments on the face of a cube for site and bond percolation, shown in
Fig. S3 of the Supplemental Information [36], we find universal behavior with
exponent $\zeta=2.25(1)$. This further confirms the vanishing cluster count
exponent in this case and provides an example of a critical system for which
$\zeta\neq 2$.

## References

  * Komogortsev  _et al._ [2019] S. V. Komogortsev, R. S. Iskhakov, and V. A. Felâk, Fractal dimension effect on the magnetization curves of exchange-coupled clusters of magnetic nanoparticles, J. Exp. Theor. Phys. 128, 754 (2019).
  * Tailleur and Cates [2008] J. Tailleur and M. E. Cates, Statistical mechanics of interacting run-and-tumble bacteria, Phys. Rev. Lett. 100, 218103 (2008).
  * Buttinoni  _et al._ [2013] I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger, and T. Speck, Dynamical clustering and phase separation in suspensions of self-propelled colloidal particles, Phys. Rev. Lett. 110, 238301 (2013).
  * Palacci  _et al._ [2013] J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine, and P. M. Chaikin, Living crystals of light-activated colloidal surfers, Science 339, 936 (2013).
  * Be’er  _et al._ [2020] A. Be’er, B. Ilkanaiv, R. Gross, D. B. Kearns, S. Heidenreich, M. Baer, and G. Ariel, A phase diagram for bacterial swarming, Commun. Phys. 3, 10.1038/s42005-020-0327-1 (2020).
  * Swart [2002] G. W. Swart, Activated leukocyte cell adhesion molecule (CD166/ALCAM): Developmental and mechanistic aspects of cell clustering and cell migration, Eur. J. Cell Biol. 81, 313 (2002).
  * Guo  _et al._ [2002] H. Guo, I. Lee, M. Kamar, S. K. Akiyama, and M. Pierce, Aberrant N-glycosylation of $\beta$1 integrin causes reduced $\alpha$5$\beta$1 integrin clustering and stimulates cell migration, Cancer Res. 62, 6837 (2002).
  * Vicsek  _et al._ [1995] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Novel type of phase transition in a system of self-driven particles, Phys. Rev. Lett. 75, 1226 (1995).
  * Cavagna and Giardina [2014] A. Cavagna and I. Giardina, Bird flocks as condensed matter, Annu. Rev. Condens. Matter Phys. 5, 183 (2014).
  * Katz  _et al._ [2011] Y. Katz, K. Tunstrøm, C. C. Ioannou, C. Huepe, and I. D. Couzin, Inferring the structure and dynamics of interactions in schooling fish, Proc. Natl. Acad. Sci. 108, 18720 (2011).
  * Gordon [2014] D. M. Gordon, The ecology of collective behavior, PLOS Biology 12, e1001805 (2014).
  * Kovács  _et al._ [2014a] I. A. Kovács, E. M. Elçi, M. Weigel, and F. Iglói, Corner contribution to cluster numbers in the Potts model., Phys. Rev. B 89, 064421 (2014a).
  * Stauffer and Aharony [1992] D. Stauffer and A. Aharony,  _Introduction to percolation theory_ (Taylor & Francis, London, 1992).
  * Smirnov [2001] S. Smirnov, Critical percolation in the plane: conformal invariance, Cardy’s formula, scaling limits, C. R. Acad. Sci. Paris Sér. I Math. 333, 239 (2001).
  * Tsai  _et al._ [2013] J. Tsai, S. C. P. Yam, and W. Zhou, Conformal invariance of the exploration path in 2-d critical bond percolation in the square lattice (2013), arXiv:1112.2017 [math.PR] .
  * Cardy [1987] J. L. Cardy, Conformal invariance, in  _Phase Transitions and Critical Phenomena_ , Vol. 11, edited by C. Domb and J. Lebowitz (Academic Press, London, 1987) p. 55.
  * Schramm [2000] O. Schramm, Scaling limits of loop-erased random walks and uniform spanning trees, Israel J. Math. 118, 221 (2000).
  * Smirnov and Werner [2001] S. Smirnov and W. Werner, Critical exponents for two-dimensional percolation, Math. Res. Lett. 8, 729 (2001).
  * Dotsenko and Fateev [1984] V. Dotsenko and V. Fateev, Conformal algebra and multipoint correlation functions in 2D statistical models, Nucl. Phys. B 240, 312 (1984).
  * Cardy [1992] J. L. Cardy, Critical percolation in finite geometries, J. Phys. A 25, L201 (1992).
  * Yu  _et al._ [2008] R. Yu, H. Saleur, and S. Haas, Entanglement entropy in the two-dimensional random transverse field ising model, Phys. Rev. B 77, 140402 (2008).
  * Kovács  _et al._ [2012] I. A. Kovács, F. Iglói, and J. Cardy, Corner contribution to percolation cluster numbers., Phys. Rev. B 86, 214203 (2012).
  * Berthiere and Witczak-Krempa [2022] C. Berthiere and W. Witczak-Krempa, Entanglement of skeletal regions, Phys. Rev. Lett. 128, 240502 (2022).
  * Kovács and Iglói [2014] I. A. Kovács and F. Iglói, Corner contribution to percolation cluster numbers in three dimensions., Phys. Rev. B 89, 174202 (2014).
  * Kovács  _et al._ [2014b] I. A. Kovács, J.-C. A. d’Auriac, and F. Iglói, Excess entropy and central charge of the two-dimensional random-bond potts model in the large-q limit., J. Stat. Mech. Theory Exp. 2014, P09019 (2014b).
  * Witczak-Krempa [2019] W. Witczak-Krempa, Entanglement susceptibilities and universal geometric entanglement entropy, Phys. Rev. B 99, 075138 (2019).
  * Cardy and Peschel [1988] J. L. Cardy and I. Peschel, Finite-size dependence of the free energy in two-dimensional critical systems, Nuclear Physics B 300, 377 (1988).
  * Cardy [1989] J. L. Cardy, Boundary conditions, fusion rules and the Verlinde formula, Nuclear Physics B 324, 581 (1989).
  * Stéphan and Dubail [2013] J.-M. Stéphan and J. Dubail, Logarithmic corrections to the free energy from sharp corners with angle 2$\pi$, J. Stat. Mech , P09002 (2013).
  * van den Berg and Conijn [2016] J. van den Berg and R. Conijn, The expected number of critical percolation clusters intersecting a line segment, Electron. Commun. Probab. 21, 1 (2016).
  * Cardy [2001] J. Cardy, Conformal invariance and percolation (2001), arXiv:math-ph/0103018 [math-ph] .
  * Kovács  _et al._ [2014c] I. A. Kovács, E. M. Elçi, M. Weigel, and F. Iglói, Corner contribution to cluster numbers in the potts model, Phys. Rev. B 89, 064421 (2014c).
  * Lorenz and Ziff [1998a] C. D. Lorenz and R. M. Ziff, Precise determination of the bond percolation thresholds and finite-size scaling corrections for the sc, fcc, and bcc lattices, Phys. Rev. E 57, 230 (1998a).
  * Lorenz and Ziff [1998b] C. D. Lorenz and R. M. Ziff, Universality of the excess number of clusters and the crossing probability function in three-dimensional percolation, J. Phys. A 31, 8147 (1998b).
  * Wang  _et al._ [2013] J. Wang, Z. Zhou, W. Zhang, T. M. Garoni, and Y. Deng, Bond and site percolation in three dimensions, Phys. Rev. E 87, 052107 (2013).
  * [36] See Supplemental Material at [LINK] for Figures S1-S4.
  * Kasteleyn and Fortuin [1969] P. W. Kasteleyn and C. M. Fortuin, Phase transitions in lattice systems with random local properties, J. Phys. Soc. Japan 26, 11 (1969).
  * Kovács and Iglói [2012a] I. A. Kovács and F. Iglói, Universal logarithmic terms in the entanglement entropy of 2d, 3d and 4d random transverse-field ising models., EPL 97, 67009 (2012a).
  * Wu [1982] F. Y. Wu, The Potts model, Rev. Mod. Phys. 54, 235 (1982).
  * Cardy [2005] J. Cardy, SLE for theoretical physicists, Ann. Phys. 318, 81 (2005), special Issue.
  * Kovács and Iglói [2012b] I. A. Kovács and F. Iglói, Universal logarithmic terms in the entanglement entropy of 2d, 3d and 4d random transverse-field ising models, EPL 97, 67009 (2012b).