Synchronization of Network-Coupled Oscillators with Uncertain Dynamics

TL;DR

This paper develops a probabilistic framework to analyze how uncertainty in natural frequencies affects synchronization in complex networks of oscillators, revealing robustness factors and network properties that promote synchronization.

Contribution

It extends the Synchrony Alignment Function (SAF) to account for uncertain, heterogeneous natural frequencies using probability theory, providing new insights into synchronization robustness.

Findings

Uncertainty in natural frequencies can enhance synchronization robustness.

Network properties like degree-frequency correlation promote synchronization.

The extended SAF framework predicts synchronization behavior under uncertainty.

Abstract

Synchronization of network-coupled dynamical units is important to a variety of natural and engineered processes including circadian rhythms, cardiac function, neural processing, and power grids. Despite this ubiquity, it remains poorly understood how complex network structures and heterogeneous local dynamics combine to either promote or inhibit synchronization. Moreover, for most real-world applications it is impossible to obtain the exact specifications of the system, and there is a lack of theory for how uncertainty affects synchronization. We address this open problem by studying the Synchrony Alignment Function (SAF), which is an objective measure for the synchronization properties of a network of heterogeneous oscillators with given natural frequencies. We extend the SAF framework to analyze network-coupled oscillators with heterogeneous natural frequencies that are drawn as a multivariate random vector. Using probability theory for quadratic forms, we obtain expressions for the expectation and variance of the SAF for given network structures. We conclude with numerical experiments that illustrate how the incorporation of uncertainty yields a more robust theoretical framework for enhancing synchronization, and we provide new perspectives for why synchronization is generically promoted by network properties including degree-frequency correlations, link directedness, and link weight delocalization.

Equations

$$\frac{d{\theta }_n}{dt}={\omega }_n+K\sum _{m=1}^N{A}_{nm}{H}_{nm}({\theta }_m-{\theta }_n),\text{for }n\in \{1,\dots ,N\}.$$

$$z=\frac{1}{N}\sum _{j=1}^N{e}^{i{\theta }_j},$$

$${\dot{\theta }}_n={\tilde{\omega }}_i-K{H}^{\prime }(0)\sum _{j=1}^N{L}_{nm}{\theta }_m,$$

$${\mathbf{\theta }}^{\ast }=\frac{L^{†}\tilde{\mathbf{\omega }}}{K{H}^{\prime }(0)}.$$

$$r\approx 1-\frac{\Vert {\mathbf{\theta }}^{\ast }{\Vert }^2}{2N}=1-\frac{{\tilde{\mathbf{\omega }}}^T{L}^{†T}{L}^{†}\tilde{\mathbf{\omega }}}{2N{K}^2{H}^{\prime 2}(0)}=1-\frac{{\tilde{\mathbf{\omega }}}^{T}U(S^{†}{)}^2{U}^T\tilde{\mathbf{\omega }}}{2N{K}^2{H}^{\prime 2}(0)},$$

$$J(\tilde{\mathbf{\omega }},L)=\frac{{\tilde{\mathbf{\omega }}}^{T}U(S^{†}{)}^2{U}^T\tilde{\mathbf{\omega }}}{N}$$

$$J(\tilde{\mathbf{\omega }},L)=\frac{1}{N}\sum _{j=2}^N\frac{({\mathbf{u}}^{jT}\tilde{\mathbf{\omega }}{)}^2}{s_j^2},$$

$$r\approx 1-\frac{J(\tilde{\mathbf{\omega }},L)}{2K^2{H}^{\prime 2}(0)},$$

$$r\approx 1-\frac{J(\mathbf{\omega },L)}{2K^2}.$$

$$\mathbb{E}\left[{\mathbf{y}}^{T}X\mathbf{y}\right]$$

$$\mathbb{E}[J(\mathbf{\omega },L)]=J(\mathbf{\mu },L)+\frac{1}{N}\mathrm{Tr}[U(S^{†}{)}^2{U}^T\Sigma ].$$

$$\mathbb{E}[J(\mathbf{\omega },L)]$$

$$\mathbb{E}[J(\mathbf{\omega },L)]=J(\mathbf{\mu },L)+\frac{\mathrm{Var}(\omega )}{N}\sum _{j=2}^N\frac{1}{s_j^2},$$

$$U(S^{†}{)}^2{U}^T=U(S^{†}{)}^T{V}^{T}V{S}^{†}{U}^T=(L^{†}{)}^T{L}^{†},$$

$$\text{Tr}[U(S^{†}{)}^2{U}^T\Sigma ]$$

$$\mathbb{E}[J(\mathbf{\omega },L)]$$

$$\mathrm{Var}\left[{\mathbf{y}}^{T}X\mathbf{y}\right]$$

$$\mathrm{Var}\left[J(\mathbf{\omega },L)\right]$$

$$0^{T}U(S^{†}{)}^2{U}^T\Sigma U(S^{†}{)}^2{U}^{T}0=0,$$

$$\mathrm{Var}\left[J(\mathbf{\omega },L)\right]$$

$$\mathrm{Tr}[[U(S^{†}{)}^2{U}^T\Sigma {]}^2]=\mathrm{Var}(\omega {)}^2\mathrm{Tr}[(L^{†}{)}^T{L}^{†}(L^{†}{)}^T{L}^{†}]=\mathrm{Var}(\omega {)}^2\sum _{j=1}^N{\hat{\lambda }}_j,$$

$${\mathbf{\mu }}^{T}U(S^{†}{)}^2{U}^T\Sigma U(S^{†}{)}^2{U}^T\mathbf{\mu }$$

$$=\mathrm{Var}(\omega ){\mathbf{\mu }}^{T}U(S^{†}{)}^4{U}^T\mathbf{\mu },$$

$$\mathrm{Var}\left[J(\mathbf{\omega },L)\right]$$

$$\ell =\frac{{\sum }_{n,m}{\chi }_{nm}{A}_{nm}}{{\sum }_{n,m}{\chi }_{nm}},$$

$${\rho }_m(\lambda )\to \rho (\lambda )=\{\begin{array}{ll}\frac{k\sqrt{4(k-1)-(\lambda -k{)}^2}}{2\pi (k^2-(\lambda -k{)}^2)},& \text{if ∣λ−k∣≤2k−1,}\\ 0,& \text{otherwise.}\end{array}$$

$${\rho }_m(\lambda )\to \rho (\lambda )=\{\begin{array}{ll}\frac{k\sqrt{4(\ell k-{\ell }^2)-(\lambda -k{)}^2}}{2\pi \ell \left[k^2-(\lambda -k{)}^2\right]},& \text{if ∣λ−k∣≤2ℓk−ℓ2,}\\ 0,& \text{otherwise,}\end{array}$$

$$\mathbb{E}[J(\mathbf{\omega },L)]$$

$$\mathrm{Var}[J(\mathbf{\omega },L)]$$

$$p_{\text{dir}}=\frac{1}{2N⟨k⟩}\sum _{n,m}\left|A_{nm}-A_{mn}\right|.$$

Full text

{AMS}

34C15, 34D06, 05C82

1 Introduction

The emergence of collective behavior in ensembles of network-coupled dynamical systems is a widely studied and active area of research in the dynamical systems and network-science communities [2, 22]. Examples where synchronization is vital for enabling robust functionality are plentiful in both natural and engineered systems. Examples of natural systems that display synchronization include cardiac pacemakers [16], circadian rhythms [46], neuronal processing [14], gene regulation [12], intestinal activity [1], and cell cycles [24]. Examples of engineered systems that require robust synchronization include Josephson junction arrays [45], electrochemical oscillators [10], and power grids [28, 31]. On the other hand, in certain systems the excess of synchronization may lead to unfavorable consequences, including bridge oscillations [41] and Parkinson’s disease, tremors, and seizures [29].

Given the widespread importance of robust synchronization in so many applications, a broad but practical question that still remains poorly understood is the following: How do the dynamical and structural properties in a network of coupled oscillators affect the macroscopic synchronization properties of the system as a whole? A large body of literature indicates that the microscopic properties of a network’s topology and the local dynamics and coupling mechanism in oscillator networks combine in very complicated ways and can lead to a wide range of novel dynamical behaviors  [9, 26, 27, 32, 33]. It is therefore natural to ask: What arrangements of network topologies and local oscillator dynamics optimize network synchronization?

In previous work, we addressed this question by developing a framework that introduced the Synchrony Alignment Function (SAF) [36, 37], which is an objective function that quantifies the interplay between complex network structure and local dynamics. The SAF quantifies the alignment of the oscillators’ heterogeneities with the network heterogeneity, as manifest in the singular vectors (or eigenvectors in the case of undirected graphs) of the combinatorial graph Laplacian matrix. In section 2.2, we present the generalized version of the SAF [37], which allows the network to possibly contain directed links. Importantly, the SAF gives rise to a nonnegative scalar that serves as an objective measure of the synchronization properties of an oscillator network and provides an avenue for optimization; it allows for the systematic optimization of synchronization properties for networks of heterogeneous oscillators under a wide range of practical constraints [36, 37, 43, 37]. For example, the SAF framework has been used to rank network links according to their importance to synchronization (which for heterogeneous oscillators, depends both on the network and the oscillator frequencies), thereby identifying network modification that judiciously enhance synchronization [43]. It has been applied to optimize heterogeneous chaotic oscillators and was shown to effectively improve the synchronization of experimental electronic circuits [34]. The SAF framework has also served useful beyond the goal of optimization—namely, it has been used to characterized a phenomenon called erosion of synchronization, which arises under frustrated coupling [38, 40].

The SAF quantifies the alignment of network heterogeneity with the oscillators’ heterogeneity—that is, their different natural frequencies $\{\omega_{i}\}$—which requires knowledge of the precise network and the particular frequencies. However, such information can be impossible to obtain for many applications. For example, any attempt to measure the oscillation frequency for a brain region, heart tissue, intestinal tissue, etc. will inevitably produce some measurement error. Understanding the effects of uncertainty on synchronization, and its optimization, remains an important open problem. Thus motivated, in this paper we relax the assumption that the frequencies are deterministic parameters. Instead, we extend the SAF framework by allowing the oscillators’ natural frequencies to be encoded by a multivariate random variable $\bm{\omega}$ with known expectation $\bm{\mu}$ and covariance matrix $\Sigma$. Using the observation that the SAF is a quadratic form of a nonnegative, symmetric matrix [25], we utilize the extensive literature on quadratic forms for random vectors to obtain expressions for the expectation and variance of the SAF. Our results provide theoretical justification for numerical experiments presented in [43] (see figure 6.2(c)), where we found the SAF to effectively optimize systems with moderate uncertainty about the frequencies. We also report new and unexpected findings, for instance that for independent and identically distributed frequencies, uncertainty always increases the SAF, in expectation, thereby inhibiting synchronization.

In the simplest case, each frequency $\omega_{i}$ is an identically and independently distributed random variable, implying $\bm{\mu}=\overline{\mu}[1,\dots,1]^{T}$ and $\Sigma=\sigma^{2}I$ where $\overline{\mu}\in\mathbb{R}$ and $\sigma^{2}\in\mathbb{R}_{+}$ are constants. This scenario describes when one has information about the network structure and the distribution of frequencies, but no knowledge about how the oscillators’ frequencies are arranged across the network. We use this setting to shed light on another important topic in the field of synchronization: What generic structural and/or dynamical properties of oscillator networks promote or inhibit synchronization? We begin by demonstrating that synchronization is promoted when link weights are delocalized, i.e., when a network has more links that are on average weaker compared to fewer links that are on average stronger. We support this finding through both numerical experiments and analytically, by combining our probabilistic analysis of the SAF with random matrix theory that predicts the spectral properties of the combinatorial graph Laplacian matrix $L$. We find that synchronization is enhanced by a concentration (i.e., lack of heterogeneity) of the spectral density for $L$ (singular values in the case of directed networks and eigenvalues in the case of undirected networks), which establishes a new connection between synchronization theory for heterogeneous oscillators and identical oscillators [3, 19, 21, 42]. In particular, the eigenratio — which is one measure for spectral heterogeneity — is well known to be important for identical oscillators [3]. Here, we find for IID frequencies that the expectation and variance of the SAF are proportional to the second and fourth moments of the spectral density, respectively, which are two different measures for spectral heterogeneity. We also show that increasing the directedness of a network (that is, the fraction of links that are directed versus undirected) promotes synchronization in a similar way as delocalization. Finally, we consider correlations between the node degrees and the oscillator frequencies, and complement previous research on the effects of these correlations [5, 20, 23, 36, 37] by exploring how different types of correlations can be used to improve or diminish synchronization.

The remainder of this paper is organized as follows. In Section 2, we present background information. In Section 3, we present our main analytical results, quantifying the expectation and variance of the SAF using a latent variable approach for natural frequency ensembles. We also demonstrate the accuracy of these results using numerical simulations. In Section 4, we apply these new results to identify generic network properties that promote synchronization, including link weight delocalization, directedness, and degree-frequency correlations. In Section 5, we conclude with a discussion.

2 Background Information and a Motivating Numerical Experiment

We first present Kuramoto’s phase oscillator model (Section 2.1), the SAF framework (Section 2.2), and a motivating experiment (Section 2.3).

2.1 The Network Kuramoto Model

Definition 2.1** (Network Kuramoto Phase Oscillator Model [11]).**

Consider a collection of $N$ phase oscillators in which $\theta_{n}\in[0,2\pi)$ is the phase of oscillator $n$, $\omega_{n}\in\mathbb{R}$ is the natural frequency of oscillator $n$, the adjacency matrix $A$ encodes the network-coupling of oscillators so that $A_{nm}$ is the weight of the link $m\to n$, $H_{nm}:\mathbb{R}\to\mathbb{R}$ is an interaction-specific, $2\pi$-periodic coupling function that is differentiable at [math], and $K>0$ is the global coupling strength. The network Kuramoto phase oscillator model [11] is given by the following system of first-order nonlinear differential equations:

[TABLE]

Kuramoto’s phase oscillator model was originally derived to describe the synchronization of weakly-coupled limit-cycle oscillators. This refers to oscillators whose steady-state behavior in the absence of coupling are stable limit cycles, and when coupling is added the geometry of the limit cycles are not significantly altered. Moreover, it is often assumed that the oscillator interactions follow a uniform functional form, i.e., $H_{nm}(\theta)=H(\theta)$ for all $n,m=1,\dots,N$. The choice $H(\theta)=\sin(\theta)$ is particularly popular and well representative of the generic case, using the first-order term of a Fourier expansion for an odd function $H(\theta)$, and is often referred to simply as the network Kuramoto model. Of particular importance is the network structure encoded by adjacency matrix $A$. The entries of $A$ indicate the existence and weight of each link in the network, such that $A_{nm}\geq 0$ is the weight of the link $m\to n$. (If no link exists, then $A_{nm}=0$.) We assume that the network is weighted, has non-negative links, and no self loops. In the case of a directed network, we define the (weighted) in- and out-degrees by $k_{n}^{\text{in}}=\sum_{m=1}^{N}A_{nm}$ and $k_{n}^{\text{out}}=\sum_{m=1}^{N}A_{mn}$, respectively. For undirected networks, one has $A_{nm}=A_{mn}$, implying $A^{T}=A$, and we define $k_{n}^{\text{in}}=k_{n}^{\text{out}}=k_{n}$. We also define the mean degree $\langle k\rangle=(1/N)\sum_{n=1}^{N}k_{n}$.

Next we define the classical Kuramoto order parameter, which is an objective measure of the degree of synchronization in the network Kuramoto model.

Definition 2.2** (The Kuramoto Order Parameter [11]).**

Consider the network Kuramoto phase oscillator model given in Definition 2.1. The Kuramoto order parameter is the complex number $z\in\mathbb{C}$ given by the centroid of all $N$ oscillator phases $\{\theta_{j}\}$ when mapped onto the complex unit circle:

[TABLE]

*where $i=\sqrt{-1}$ is the imaginary unit. *

It is convenient to use the polar decomposition of the order parameter $z=re^{i\psi}$, where $r=|z|\in[0,1]$ is the magnitude of $z$ and $\psi$ is the collective phase of the system. In particular, $r$ is an objective measure of the degree of synchronization since an incoherent state yields $r\approx 0$, and a strongly synchronized state yields $r\approx 1$.

2.2 The Synchrony Alignment Function

Next, we summarize the SAF framework, which was first introduced in [36] for undirected networks and later generalized for directed networks [37]. To help facilitate the presentation of our main analytical results in section 3, we will define the SAF in terms of a quadratic form [25]. The fact that the SAF takes this form has not been noted in previous research.

We begin with the assumption that the coupling function is uniform, $H_{nm}(\theta)=H(\theta)$, and that $H(\theta)$ is smooth and positively-sloped at zero, and $H(0)\ll 1$ with a root $H(\theta)=0$ nearby. As shown/discussed in [38], these assumptions allow for the possibility of a strongly synchronized state when the coupling strength $K$ is for a sufficiently large or when the spread of oscillator frequencies is sufficiently small. We then consider the strongly synchronized regime where phase angles satisfy $|\theta_{m}-\theta_{n}|\ll 1$, allowing us to linearize Eq. 1 to obtain the linearized system

[TABLE]

where $\tilde{\omega}_{n}=\omega_{n}+KH(0)k_{n}^{in}$ is the effective natural frequency of oscillator $n$ and $L$ is the combinatorial Laplacian matrix whose entries are given by $L_{nm}=\delta_{nm}k_{n}^{in}-A_{nm}$.

The combinatorial graph Laplacian $L$ is a cornerstone for numerous problems in network science, and we will make use of its singular value decomposition. The singular value decomposition of $L$ is given by $L=USV^{T}$, where $S=\text{diag}(s_{1},s_{2},\dots,s_{N})$, and the columns of $U$ and $V$ are populated by the orthonormal left and right singular vectors $\{\bm{u}^{i}\}_{i=1}^{N}$ and $\{\bm{v}^{i}\}_{i=1}^{N}$, which satisfy $L\bm{v}^{i}=s_{i}\bm{u}^{i}$ and $L^{T}\bm{u}^{i}=s_{i}\bm{v}^{i}$ [8]. For undirected networks $L$ is symmetric, $L=L^{T}$, and the singular values and singular vectors are identical to the eigenvalues and eigenvectors of $L$, and so it suffices to consider only the eigenvalue decomposition is required $L=V\Lambda V^{T}$, where $\Lambda=\text{diag}(\lambda_{1},\lambda_{2},\dots,\lambda_{N})$ and the columns of $V$ are populated by the orthonormal eigenvectors.

The SAF framework utilizes the pseudo-inverse of $L$, which is generally defined by $L^{\dagger}=VS^{\dagger}U^{T}$, where $S^{\dagger}=\text{diag}(0,s_{2}^{-1},\dots,s_{N}^{-1})$ [4]. For undirected networks one has $L^{\dagger}=V\Lambda^{\dagger}V^{T}$ with $\Lambda^{\dagger}=\text{diag}(0,\lambda_{2}^{-1},\dots,\lambda_{N}^{-1})$. Note that we explicitly assume there are $N-1$ singular values and nonzero eigenvalues, which occurs, for example, if the network is strongly connected. Under this assumption, the real part of all other eigenvalues is positive, $\text{Re}(\lambda_{i})>0$ for $i=2,\dots,N$ and all other singular values are positive and can be ordered $0=s_{1}<s_{2}\leq\cdots\leq s_{N}$.

Returning to the dynamics of Eq. 3, we search for a phase-locked solution of the form $\theta_{n}(t)=\theta_{n}^{*}+\Omega t$, where $\Omega$ is the collective frequency of the synchronized state. In the case of an undirected network we have $\Omega=\langle\tilde{\bm{\omega}}\rangle$, whereas in the directed case we have $\Omega=(\bm{u}^{1T}\tilde{\omega})/(\bm{u}^{1T}\bm{1})$ [39]. In either case, we enter the rotating reference frame $\theta_{n}\mapsto\theta_{n}+\Omega t$, and search for the phase-locked state in this rotating frame satisfying $\dot{\theta}_{n}=0$ for all $n=1,\dots,N$. The solution is given, in vector form, by

[TABLE]

Returning to the Kuramoto order parameter in Eq. 2, in the strongly synchronized regime we have that

[TABLE]

where we have used that $V^{T}V=I$ and $S^{\dagger}$ is diagonal so that $(S^{\dagger})^{T}S^{\dagger}=(S^{\dagger})^{2}=\text{diag}(0,s_{2}^{-2},\dots,s_{N}^{-2})$. Equation Eq. 5 brings us to the quadratic form definition of the SAF

Definition 2.3** (Synchrony Alignment Function (SAF) [36, 37]).**

Consider the network Kuramoto phase oscillator model given in Definition 2.1 with linearized dynamics given by Eq. 3 with effective natural frequencies $\tilde{\bm{\omega}}$ and Laplacian $L$. The Synchrony Alignment Function is a function $J:\mathbb{R}^{N}\times\mathbb{M}_{N\times N}\to\mathbb{R}$ defined as

[TABLE]

*where $(S^{\dagger})^{2}=\text{diag}(0,s_{2}^{-2},\dots,s_{N}^{-2})$ is populated by the inverse squares of the nontrivial singular values of $L$ and the columns of $U$ are populated by the left singular vectors $\{\bm{u}^{i}\}_{i=1}^{N}$ of $L$. *

For the remainder of this paper, we will refer to (6) simply as the SAF. We emphasize that the SAF is given here as a quadratic form with matrix $U(S^{\dagger})^{2}U^{T}$. It appears different from, although is equivalent to, the forms given in previous work [36, 37]. In particular, one can alternately write the pseudo-inverse of the Laplacian as $L^{\dagger}=\sum_{j=2}^{N}\bm{v}^{j}\bm{u}^{jT}/s_{j}$, which leads to the alternate form

[TABLE]

which more explicitly describes the SAF as a weighted combination of projections of the frequency vector $\tilde{\bm{\omega}}$ onto the nontrivial left singular vectors $\bm{u}^{j}$. However, the quadratic form in Eq. 6 has some advantages that we will leverage below in our analysis of the expected value and variance of $J(\tilde{\bm{\omega}},L)$ when considering randomized frequencies.

Together, Eq. 5 and Eq. 6 yield

[TABLE]

thus allowing us to use the SAF as an objective measure for the degree of synchronization in a network. On its own, the SAF quantifies the interplay between a network’s dynamics (i.e., through the natural frequencies $\{\omega_{i}\}$) and its structure (i.e., through Laplacian singular values $\{s_{i}\}$ and singular vectors $\{\bm{u}_{i}\}$ and $\{\bm{v}_{i}\}$). Equation Eq. 8 implies that, in the linearized regime where the approximations made to obtain Eq. 3 are valid, the synchronization properties of a given network are promoted when the SAF $J(\tilde{\bm{\omega}},L)$ is made small. Specifically, by inspecting Eq. 7, we can see that this is done by aligning the frequency vector with the singular vectors associated with large singular values and making it orthogonal to the singular vectors associated with small singular values. We also note here that in the typical case of sinusoidal coupling, $H(\theta)=\sin(\theta)$, we have that $H^{\prime}(0)=1$ and $H(0)=0$, implying $\tilde{\omega}_{n}=\omega_{n}$ and simplifying Eq. 8 to

[TABLE]

In the remaining text, we will assume sinusoidal coupling, $H(\theta)=\sin(\theta)$.

2.3 Motivating Numerical Experiment

The SAF allows us to evaluate the synchronization properties of a specific network with a specific set of natural frequencies. Moreover, Eq. 8 allows for the efficient optimization of a network’s synchronization properties under various constraints. As an illustrative example, we consider a directed Erdős-Rényi (ER) network [7] of size $N=500$ nodes with a mean degree $\langle k\rangle=10$ and a given set of $N$ natural frequencies $\{\omega_{i}\}$ drawn identically and independently from a normal distribution with zero mean and unit variance. We consider several different arrangements (i.e., permutations) of the given set of frequencies, a range of coupling strengths $K\in[0,0.5]$ and. For each arrangement and $K$ value, we numerically integrate Eq. 1 with sinusoidal coupling and a random initial condition using Newton’s forward-difference method.

In Fig. 1 (a), we plot the degree of synchronization given by the order parameter $r$ as a function of the coupling strength $K$ for $40$ different arrangements of the same collection of natural frequencies on the network. To demonstrate the utility of the SAF we find the minimum and maximum value of $J(\bm{\omega},L)$ over these $40$ arrangements (given by $0.0158$ and $0.0228$, respectively), and we plot the corresponding curves in blue and red, respectively. As predicted by Eq. 8 and Eq. 9, the minimum and maximum values of the SAF correspond to the maximum and minimum values of the order parameter, respectively, in the strongly synchronized regime. This is highlighted in the inset where we zoom-in on the synchronized regime. Note that once synchronization occurs, the precise value of the order parameter varies significantly from one arrangement to the next. To highlight this, in Fig. 1 (b) we plot an empirically observed distribution of $r$ — which we denote $f(r)$ — that is estimated using 5000 random arrangements at $K=0.25$, i.e., near the beginning of the strongly synchronized regime. Despite the fact that these results come from a single network structure and a single set of natural frequencies, the distribution of observed order parameter values is remarkably wide, with a bulk ranging from roughly $0.74$ to $0.82$. (The sample mean and standard deviation are $0.7861$ and $0.0156$, respectively.)

As illustrated with this example, the SAF framework has proven very useful in both characterizing and optimizing the synchronization properties of coupled oscillator networks with specific information about the network structure and local dynamics, i.e., natural frequencies $\{\omega_{i}\}$. However, the requirement knowing the precise natural frequencies introduces a hurdle in quantifying and developing an understanding of the overall synchronization properties of, for instance, a network structure without specific knowledge of the natural frequencies. Instead, some statistical properties of the natural frequencies in a network are more likely attainable, in which case a description of the expected value of the SAF and the likely size of any deviations, i.e., its variance, is a more realistic outcome. Moreover, this approach lends itself better to providing a broader understanding of the network properties that promote or inhibit synchronization. To this end, in the following section we will develop theory for the expectation and variance of the SAF by treating the natural frequencies as latent random variables.

3 Analysis of the SAF under Uncertain Frequencies

We will now study the behavior of the SAF in a new context where natural frequencies are not determined, but rather are random variables. We consider ensembles of natural frequency vectors and describe the expectation (Section 3.1) and variance (Section 3.2) of the SAF. In (Section 3.3), we provide numerical support to validate these results.

Before continuing, we emphasize that we will still treat the network structure as fixed, and therefore the quantities that result from the network structure (e.g., the Laplacian $L$, singular vectors $\bm{u}^{i}$ and $\bm{v}^{i}$, singular values $s_{i}$, etc.) are also fixed. On the other hand, from this point forward, we will treat each natural frequency $\omega_{i}$ as a random variable with mean $\mu_{i}$ and variance $\sigma_{i}^{2}$. The vector $\bm{\omega}=[\omega_{1},\dots,\omega_{N}]^{T}$ is then a vector of random variables with mean vector $\bm{\mu}=[\mu_{1},\dots,\mu_{N}]^{T}$ and covariance matrix $\Sigma$ (whose diagonal entries are the variances $\sigma_{i}^{2}$ and whose off-diagonal entries are the covariances). Therefore, $\bm{\mu}$ and $\Sigma$ may be interpreted as latent variables for the natural frequencies $\bm{\omega}$.

3.1 Expectation of the SAF under Uncertain Frequencies

We first study the expectation $\mathbb{E}[J(\bm{\omega},L)]$ of the SAF, and we will present several results for different assumptions about $\bm{\mu}$ and $\Sigma$. Our analysis leverages the observation that the SAF given in Definition 2.3 is a quadratic form, allowing us to utilize known results for the quadratic forms of random vectors.

Theorem 3.1** (Expectation for Quadratic Forms of Random Vectors [25]).**

Consider a random vector $\bm{y}=[y_{1},\dots,y_{N}]^{T}$ in which $\bm{\mu}=[\mu_{1},\dots,\mu_{N}]^{T}$ encodes the entries’ means and $\Sigma$ encodes the covariances between entries. For any real-valued square, symmetric matrix $X$, the mean of the quadratic form $\bm{y}^{T}X\bm{y}$ is given by

[TABLE]

*where $\mathrm{Tr}[\cdot]$ denotes the matrix trace. *

Proof 3.2**.**

*This result is well known and is given as Theorem 5.2a in [25]. To provide the reader insight, we include a short proof in Appendix A. *

We now utilize Theorem 3.1 to find the expectation of the SAF for a number of different choices of random frequency vectors. First we consider the SAF for general frequency vectors $\bm{\omega}$ with given mean vector $\bm{\mu}$ and covariance matrix $\Sigma$.

Theorem 3.3** (Expected SAF with Random Frequencies).**

Consider the network Kuramoto phase oscillator model given in (2.1) with linearized dynamics given by Eq. 4 for a Laplacian matrix $L$ and natural frequencies that are drawn as a random vector $\bm{\omega}$ in which the entries’ means are given by $\bm{\mu}=[\mu_{1},\dots,\mu_{N}]^{T}$ and their covariances are given by $\Sigma$. The expected value of the SAF is given by

[TABLE]

Proof 3.4**.**

*Considering the quadratic form in Eq. 6, we let $X=N^{-1}U(S^{\dagger})^{2}U^{T}$ and $\bm{y}=\bm{\omega}$, which, when substituted into Eq. 10, simplifies to Eq. 11. *

The general form of the expected value of the SAF presented in Theorem 3.3 deserves a number of remarks. Overall, the expectation of the SAF under random frequencies can be viewed as a perturbation of the SAF computed using the expected frequencies, $J(\bm{\mu},L)$. Thus, in the limit as the variance of each natural frequency approaches zero, we recover $\mathbb{E}[J(\bm{\omega},L)]=J(\bm{\mu},L)$. Next we present several corollaries of Theorem 3.3 that consider the specific cases of natural frequency vector with equal means, independent natural frequencies with equal variances, and finally independently and identically distributed (IID) natural frequencies.

Corollary 3.5** (Expected SAF for Zero-Mean Random Frequencies).**

In the case where the natural frequencies have zero mean, $\bm{\mu}=\bm{0}$, then the expected value of the SAF is given by

[TABLE]

Proof 3.6**.**

*The result is immediate using $J({\bf 0},L)=0$. *

Corollary 3.7** (Expected SAF for Independent Random Frequencies of Equal Variance).**

In the case where the natural frequencies are uncorrelated and have the same variance, $\Sigma=\mathrm{Var}(\omega)I_{N}$, where $I_{N}$ is an identity matrix of size $N$, then the expected value of the SAF is given by

[TABLE]

*where $s_{j}$ is the $j^{\text{th}}$ singular value of the Laplacian $L$. *

Proof 3.8**.**

Using the definitions of $\Sigma$ and trace, and noting that

[TABLE]

we have that

[TABLE]

*where we have used that the trace of a matrix equals the sum of its eigenvalues and we denote the eigenvalues of $(L^{\dagger})^{T}L^{\dagger}$ as $\{\hat{\lambda}_{i}\}_{i=1}^{N}$. Since $L^{\dagger}$ is the pseudo-inverse of the Laplacian $L$ we have that $\hat{\lambda}_{1}=0$ and $\hat{\lambda}_{i}=s_{i}^{-2}$ for $i=2,\dots,N$, which completes the proof. *

Corollary 3.9** (Expected SAF for IID Frequencies).**

In the case where the natural frequencies are identically and independently distributed with zero mean, then the expected value of the SAF is given by

[TABLE]

Proof 3.10**.**

*The result follows directly from applying Corollary 3.5 and Corollary 3.7. *

Note that because $\mathrm{Var}(\omega)>0$ and $\sum_{j=2}^{N}{s_{j}^{-2}}>0$, if they are IID then the uncertainty for the natural frequencies, always increases the SAF in expectation, thereby inhibiting synchronization.

3.2 Variance of the SAF for Normally Distributed Frequencies

We now study the variance of the SAF, $\mathrm{Var}[J(\bm{\omega},L)]=\mathbb{E}[J^{2}(\bm{\omega},L)]-\mathbb{E}[J(\bm{\omega},L)]^{2}$. Similar to the previous subsection, we will present several results that are obtained under different assumptions about $\bm{\mu}$ and $\Sigma$. We will again leverage the quadratic form of the SAF given in Definition 2.3, now with a theorem describing the variance of a quadratic forms of random vectors. In order to obtain analytical expressions, we now restrict our attention to the case of normally distributed natural frequencies.

Theorem 3.11** (Variance for Quadratic Forms of Normal-Distributed Random Vectors [25]).**

Consider the random vector $\bm{y}$ and matrix $X$ from Thm. 3.1, and further assume $\bm{y}\sim\mathcal{N}(\bm{\mu},\Sigma)$ is drawn from a multivariate Gaussian distribution. Then the variance of the quadratic form $\bm{y}^{T}X\bm{y}$ is given by

[TABLE]

Proof 3.12**.**

*See Theorem 5.2c of [25] for a derivation the uses the moment generating function for $\bm{y}^{T}X\bm{y}$. *

We will now utilize Theorem 3.11 to evaluate the variance of the SAF for a number of different choices of random frequency vectors. First we consider the SAF for general frequency vectors $\bm{\omega}$ with given mean vector $\bm{\mu}$ and covariance matrix $\Sigma$.

Theorem 3.13** (Variance of SAF for Normal-Distributed Frequencies).**

Consider the network of oscillators with random frequencies described in Theorem 3.3 and further assume the frequency vector $\bm{\omega}\sim\mathcal{N}(\bm{\mu},\Sigma)$ is drawn from a multivariate Gaussian distribution with mean $\bm{\mu}$ and variance $\Sigma$. Then the variance of the SAF $J(\bm{\omega},L)$ is given by

[TABLE]

Proof 3.14**.**

*Similarly as the proof to Theorem 3.3, we use the SAF given via a quadratic form in Eq. 6, and let $X=N^{-1}U(S^{\dagger})^{2}U^{T}$ and $\bm{y}=\bm{\omega}$. Applying Theorem 3.13 yields the desired result. *

The variance of the SAF presented in Theorem 3.13 does not have as straight forward of an interpretation as the expected value (see Theorem 3.3). However, it has an important similarity in that it consists of two terms, one term is a quadratic form with the mean vector $\bm{\mu}$, and the other term is a trace of a matrix product that includes the covariance matrix $\Sigma$. Next we present several corollaries of Theorem 3.3 that, similar to the results given in section 3.1, consider the specific cases of natural frequency vector with equal means, independent natural frequencies with equal variances, and IID natural frequencies.

Corollary 3.15** (Variance of SAF for Equal-Mean Random Frequencies).**

In the case where the natural frequencies have the same mean, $\bm{\mu}=\overline{\omega}{\bf 1}$, then the variance of the SAF is given by

[TABLE]

Proof 3.16**.**

Similarly to the proof for Corollary 3.5, we may rescale the mean vector to $\bm{\mu}=\bm{0}$. The latter term on the right hand side of Eq. 18 then vanishes, i.e.,

[TABLE]

*completing the proof *

Corollary 3.17** (Variance of SAF for Independent Random Frequencies of Equal Variance).**

In the case where the natural frequencies are uncorrelated and have the same variance, $\Sigma=\mathrm{Var}(\omega)I_{N}$, where $I_{N}$ is an identity matrix of size $N$, then the variance of the SAF is given by

[TABLE]

Proof 3.18**.**

Applying Theorem 3.13 and beginning with the first term on the right hand side of Eq. 18, we have that

[TABLE]

where we have used that $U(S^{\dagger})^{2}U^{T}=(L^{\dagger})^{2}L^{\dagger}$ (see Eq. 14) and $\mathrm{Tr}[(L^{\dagger})^{T}L^{\dagger}(L^{\dagger})^{T}L^{\dagger}]$ is equal to the sum of its eigenvalues, denoted $\{\hat{\lambda}_{j}\}_{j=1}^{N}$. Moreover, similar to the proof for Corollary 3.7 we have that $\hat{\lambda}_{1}=0$ and $\hat{\lambda}_{i}=s_{i}^{-4}$ for $i=2,\dots,N$, completing the first term on the right hand side of Eq. 21. For the other term, we take the second term on the right hand side of Eq. 18, which simplifies to

[TABLE]

*completing the second term on the right hand side of Eq. 21, completing the proof. *

Corollary 3.19** (Variance of SAF for IID Frequencies).**

In the case where the natural frequencies identically and independently distributed, then the variance of the SAF is given by

[TABLE]

Proof 3.20**.**

*The result follows directly from applying Corollary 3.15 and Corollary 3.17. *

3.3 Numerical Validation of Theory and Definition of Link Weight Localization

We conclude this section by presenting numerical experiments to illustrate the accuracy of the analytical results developed above. In an attempt to kill two birds with one stone, we will present an experiment that both validates our theory and introduces a network property that will be explored in detail in the next section: link weight localization. Roughly speaking, weight localization describes the trade-off between sparser network structures with more strongly weighted links versus denser network structures with more weakly weighted links. A reasonable metric for measuring weight localization of a given network is the mean weight associated with each link in the network.

Definition 3.21** (Weight Localization).**

Consider a weighted network consisting of $N$ nodes with non-negative adjacency matrix $A$. The weight localization, denoted $\ell$ is given by the mean link weight,

[TABLE]

*where $\chi$ is the link indicator function, i.e., $\chi_{nm}=1$ if $A_{nm}>0$ and [math] if $A_{nm}=0$. *

We will use link weight localization to define an ensemble of Kuramoto phase oscillator systems, allowing us to illustrate the accuracy of our theory for a range of SAF values. We consider weighted versions of Erdős–Rényi (ER) and scale-free (SF) networks constructed using the configuration model [17]. Specifically, given a network with $N$ nodes we construct and $M=N\langle k\rangle/\ell$ undirected edges, each having weight $\ell$. (Recall that $k_{i}=\sum_{j}A_{ij}$ denotes the strength, or weighted degree, of node $i$ and $\langle k\rangle$ denotes the average of $k_{i}$.)

We use this class of networks to illustrate the utility of the SAF in predicting the Kuramoto order parameter $r$ in the strong synchronization regime and verify the analytical results from the previous section. We consider $100$ weighted ER and SF networks each, all of size $N=500$ with mean degree $\langle k\rangle=10$ (and $\gamma=3$ for the SF networks). For each network we chose $\ell$ randomly and uniformly from the interval $[0.1,2]$ and consider $50$ different frequency vectors with IID normal-distributed entries with unit variance and zero mean. For this choice of frequencies, the expectation and variance of the SAF are described by Corollary 3.9 and Corollary 3.19, respectively. In Fig. 2, we plot the time average of $r$ for $K=0.5$ versus the SAF $J(\bm{\omega},L)$ for each frequency vector and each network. The dashed line indicates the linear approximation $r\approx 1-J(\bm{\omega},L)/(2K^{2})$ given by Eq. 9. Each data point is colored according to $\ell$. For both ER and SF networks, we note strong agreement between the predicted and observed values of $r$.

Next, we verify our analysis for the expectation and variance of the SAF by comparing these predictions to sample means and sample variances that are estimated using the $50$ frequency vectors for each network. In Fig. 3 we plot these sample means and variances versus their expected SAF and its variances, which are given by Eq. 16 and Eq. 24, respectively. Again, we color each datapoint according to $\ell$. Overall, we note excellent agreement between our analytical predictions and the sample means and variances observed. We note a small deviation in the variance for significantly delocalized network, where the analytical results over-predict our observations.

4 SAF Expectation and Variance Reveals Generic Network Properties that Promote Synchronization

In this section, we illustrate that the theoretical results obtained in the previous section provide a framework for understanding how generic network properties affect synchronization. In particular, we will study link weight delocalization (Section 4.1), link directedness (Section 4.2), and degree-frequency correlations (Section 4.3), and we show that each of these properties tends to promote the synchronization properties of network-coupled heterogeneous oscillators.

4.1 Weight Delocalization Promotes Synchronization

The first network property that we identify as a promoter of network synchronization is weight delocalization. Recall from Section 3.3 that weight localization/delocalization quantifies the trade-off between sparse/dense structures and strong/weak connections. Specifically, in the context of conserving the total of all link weights in the network, one may strengthen the average link weight in the network at the expense of ending up with fewer links (and therefore a sparser structure) or one may increase the number of total links in the network (therefore ending up with a denser structure) at the expense of weakening the average link weight. At first glance the effect that this trade-off has on overall synchronization dynamics is not clear, however we will show that weight delocalization improves synchronization while weight localization diminishes it.

We begin with some numerical simulations by plotting in Fig. 4 the Kuramoto order parameter $r$ versus coupling strength $K$ for networks with different levels of localization ranging from $\ell=2$ (localized; red, bottom) to $\ell=0.25$ (delocalized; blue, top). Similar to the simulation in the previous section, we consider IID normal-distributed natural frequencies $\{\omega_{i}\}$ and (a) ER and (b) SF networks of size $N=500$ with mean degree $\langle k\rangle=10$ with exponent $\gamma=3$ for SF networks. Each datapoint represents an average of $r$ over $50$ networks. Importantly, note that in both ER and SF networks the delocalized networks (i.e., small $\ell$) display stronger synchronization properties than the localized networks (i.e., larger $\ell$). That is, synchronization is promoted by having more links with weaker link weights.

To explain this phenomenon, we first apply and extend some results from random matrix theory to understand the effect that weight delocalization has on the eigenvalue spectrum of the network’s Laplacian matrix. Note first that the results in Fig. 4 used IID natural frequencies, and therefore the analytical description of the expected value and variance of the SAF were described by Corollary 3.9 and Corollary 3.19. Moreover, in the context of undirected networks where the Laplacian $L$ is symmetric, i.e., $L^{T}=L$, the singular value decomposition reduces to the eigenvalue/eigenvector diagonalization of $L$. In particular, the singular values $\{s_{i}\}_{i=1}^{N}$ are given precisely by the eigenvalues $\{\lambda_{i}\}_{i=1}^{N}$. Therefore, to understanding the effect of weight delocalization on the expectations of the SAF, and thereby the overall synchronization properties of a network, we must understand the effect of weight delocalization on the spectral density $\rho(\lambda)$ of Laplacian eigenvalues of a network. Note that once the behavior of $\rho(\lambda)$ is well understood, both $\mathbb{E}[J(\bm{\omega},L)]$ and $\mathrm{Var}[J(\bm{\omega},L)]$ can be easily predicted since they are proportional to the variance of frequencies—which we denote $\mathrm{Var}(\omega)$—as well as the second and fourth moments of $\rho(\lambda)$ given by $\sum_{j=2}^{N}\lambda_{j}^{-2}$ and $\sum_{j=2}^{N}\lambda_{j}^{-4}$, respectively. [See Eq. 16 and Eq. 24.]

We first study the effect of localization on the quantities $\sum_{j=2}^{N}\lambda_{j}^{-2}$ and $\sum_{j=2}^{N}\lambda_{j}^{-4}$ for $k$-regular networks, which have the property that all nodes have precisely the same degree, $k_{n}=k$. Moreover, we begin with unweighted case, where each of a node’s $k$ links have weight $w=1$. In this case McKay’s theorem characterizes the eigenvalue spectrum of the Laplacian matrix for the asymptotic, large-$N$ limit.

Theorem 4.1** (McKay’s Theorem [13]).**

Consider a sequence of unweighted $k$-regular networks with Laplacian matrices $L^{(m)}$, each of size $N_{m}$ with $N_{m}\to\infty$ as $m\to\infty$ . Further, let $\{\lambda_{i}^{(m)}\}_{i=1}^{N_{m}}$ denote the $N_{m}$ eigenvalues for each $L^{(n)}$. Then the empirical spectral density $\rho_{m}(\lambda)=\frac{1}{N_{m}}\sum_{i=1}^{N_{m}}\delta_{\lambda_{i}}$ converges in probability with $m\to\infty$ as

[TABLE]

McKay’s Theorem is a pioneering result in random matrix theory and research into the behavior of spectral densities for various matrices. The density $\rho(\lambda)$, which is centered at the node degree $k=\int_{\lambda}\lambda\hat{\rho}(\lambda)d\lambda$ and has bounded support on the interval $\lambda\in(k-2\sqrt{k-1},k+2\sqrt{k-1})$. We note that there remains an active research community working to extend McKay’s law. In particular, it was recently shown that the same result can be obtained by allowing $k$ to increase with $n$, provided that $k$ grows sufficiently slowly [6, 44], in which case $\rho(\lambda)$ converges to the famous “semi-circle distribution” [15].

With Theorem 4.1 in hand, we turn our attention to weighted $k$-regular networks. Given a weight $\ell$ a $k$-regular network consists of $N$ nodes that each have $k/\ell$ links with weight $\ell$. (Note that we must choose $\ell$ so that $k/\ell$ is a positive integer less than $N$.) We now a new theorem that is akin to McKay’s Thereom but is modified for weighted $k$-regular networks.

Theorem 4.2** (Spectral Density of Weighted $k$-Regular Networks).**

Consider a sequence of weighted $k$-regular networks with weight $\ell$ with Laplacian matrices $L^{(m)}$, each of size $N_{m}$ with $N_{m}\to\infty$ as $m\to\infty$. Further, let $\{\lambda_{i}^{(m)}\}_{i=1}^{N_{m}}$ denote the $N_{m}$ eigenvalues for each $L^{(n)}$. Then the empirical spectral density $\rho_{m}(\lambda)=\frac{1}{N_{m}}\sum_{i=1}^{N_{m}}\delta_{\lambda_{i}}$ converges in probability with $m\to\infty$ as

[TABLE]

Proof 4.3**.**

See appendix B.

A key characteristic of the spectral density $\rho(\lambda)$ is its support and heterogeneity. Specifically, its support is the interval $(k-2\sqrt{\ell k-\ell^{2}},k+2\sqrt{\ell k-\ell^{2}})$. Thus, as the weight $\ell$ increases the support widens, yielding more heterogeneous eigenvalues, and as $\ell$ decreases the support decreases, yielding more homogeneous eigenvalues (i.e., spectral concentration). In Fig. 5, we plot observed spectral densities for several choices of $\ell\in\{0.2,0.3,0.5,1,2\}$ (colored blue to red) using $100$ networks of size $N=500$ with degree $k=10$ for each value of the weight. Symbols represent the overall observations and theoretical approximations from Eq. 27 are plotted as solid curves. We note an excellent agreement between them.

Importantly, given an approximation of the spectral density, we can approximate the expected value and the variance of the SAF, given in Eq. 16 and Eq. 24. For sufficiently large $N$ the expected value of the SAF can be approximated by

[TABLE]

On the other hand, the variance of the SAF can be approximated by

[TABLE]

In Fig. 6 (a) and (b), we plot the expected value and variance, respectively, of the SAF versus localization $\ell$ for each of the $100$ networks used to create Fig. 5. Observations from individual networks are plotted in circles and the dashed curves represent the predictions given by Eq. 28 and Eq. 29, respectively

These results deserve a few remarks. First, we are ultimately interested in the behavior of the expected value and variance of the SAF via the sums $\sum_{j\geq 2}\lambda_{j}^{-2}$ and $\sum_{j\geq 2}\lambda_{j}^{-4}$, respectively. Second, we note that as network weights become more delocalized the spectral density becomes more homogeneous (i.e., concentrated), with eigenvalues clustering more tightly around the mean degree $k$. This yields smaller values of $\mathbb{E}[J(\bm{\omega},L)]$ and $\mathrm{Var}[J(\bm{\omega},L)]$, i.e., better and less variable synchronization properties of a network. Also note that, because the degree $k$ remains constant, even as localization is varied, the mean of the eigenvalue distribution remains constant since the sum of the eigenvalues equals the trace of the Laplacian, which is in turn equal to the sum of the degrees, $Nk$. Thus, as the network links become more delocalized the smaller eigenvalues increase away from zero, contributing smaller values to the sums $\sum_{j\geq 2}\lambda_{j}^{-2}$ and $\sum_{j\geq 2}\lambda_{j}^{-4}$, thereby improving synchronization. On the other hand, as the network links become more localized the smaller eigenvalues move closer to zero, contributing larger values to the sums $\sum_{j\geq 2}\lambda_{j}^{-2}$ and $\sum_{j\geq 2}\lambda_{j}^{-4}$, thereby increasing the SAF and inhibiting synchronization.

We now return to ER and SF networks, showing that a similar phenomenon occurs whereby delocalization yields a concentration of the spectral density, which ultimately decreases the SAF. We illustrate this in Fig. 7 (a) and (b), where we plot the spectral densities observed from collections of $100$ ER and SF networks, respectively, for several localization values ranging from $\ell=0.2$ (blue circles) to $2$ (red squares). Even for the case of ER networks, which are weakly heterogeneous compared to SF networks, it is clear that the spectral densities differ significantly from those of $k$-regular networks plotted in Fig. 5. However, the same qualitative behavior persists in the context of varying weight localization as with the case of the $k$-regular networks. In particular, weight delocalization (localization) results in spectral densities that are more (less) clustered about the mean, which remains constant. Therefore, weight delocalization increases the smallest eigenvalues, resulting in both smaller expected values and variances of the SAF, thereby promoting synchronization.

4.2 Network Directedness Promotes Synchronization

In our analysis above we have shown that weight delocalization tends to improve the synchronization properties of networks, i.e., networks with more, weaker connections synchronize heterogeneous oscillators better than networks with fewer, stronger connections. Another fair interpretation of this phenomenon is that, to improve a network’s synchronization properties, one should seek to connect more distinct pairs of oscillators, even at the expense of making the connections weaker. We now turn to another network property that can be used to attain a similar effect, namely network directedness, which we define as follows.

Definition 4.4** (Network Directedness).**

Consider a weighted network consisting of $N$ nodes with non-negative adjacency matrix $A$. The directedness of the network, denoted $p_{\text{dir}}$ is given by fraction of weighted links $m\to n$ that do not have an equal and opposite link $n\to m$,

[TABLE]

The directedness $p_{\text{dir}}$ can also be interpreted as the normalized sum over all entries of the absolute value of the matrix $A-A^{T}$. Thus, in the case of an undirected network we have that $A-A^{T}$ is identically zero, yielding $p_{\text{dir}}=0$. On the other hand, when each and every link in the network has no opposite counterpart, i.e., each non-zero entry of $A$ corresponds to to a zero entry of $A^{T}$, we have that $p_{\text{dir}}=1$ after accounting for the normalization by $2N\langle k\rangle$. Intermediate scenarios, where some links do and some links do not have an opposite counterpart, intermediate values of $p_{\text{dir}}$ are taken.

In the context of studying network structures having the same size and degree (i.e., the same total link weight $\sum_{ij}A_{ij}$), increasing the directedness of a network increases the total number of pairs of nodes in the network that have a connection between them, similar to the process of delocalization. This is illustrated in Fig. 8, where four nodes with just one undirected link with weight $\ell=1$ are shown in panel (a). In panel (b), we illustrate a delocalized alternative, where the set of four nodes instead share two undirected links, now with with weight $\ell=1/2$. Note that this coincides with two, not one, distinct pairs of connected nodes, while the total number of weighted links is the same as in panel (a). In panel (c), we illustrate the directed alternative, where the four nodes now share two binary directed links, again yielding two distinct pairs of connected nodes as in panel (b). We note also that in certain practical scenarios networks structures may be required to be binary, making delocalization unfeasible. Tuning directedness, however provides a more realistic alternative given the ability to connect more pairs of nodes with unweighted directed links.

We proceed by investigating the effect of directedness on synchronization by plotting in Fig. 9 the Kuramoto order parameter $r$ from a collection of unweighted ER networks of size $N=500$ with mean degree $\langle k\rangle=10$ that are iteratively rewired to become more directed. In particular, we vary directedness from $p_{\text{dir}}=0$ (totally undirected, red) to $1$ (completely directed, blue). Each data point represents the average of $r$ over a collection of $50$ networks. To emphasize the improved synchronization properties of more directed networks, we also plot in the inset a zoomed-in view of the synchronization profile in the strongly synchronized regime. To better understand this behavior, we consider a collection of $500$ ER networks (also with $N=500$ and $\langle k\rangle=10$), each assigned a directedness $p_{\text{dir}}$ randomly and uniformly drawn from the interval $[0,1]$. In Fig. 10 (a) and (b), we plot the expected value and variance of the SAF, respectively. For each case we also plot the least squares line, which for the data collected here are given by $\mathbb{E}[J(\bm{\omega},L)]\approx-4.70\times 10^{-3}p_{\text{dir}}+2.11\times 10^{-2}$ and $\mathrm{Var}[J(\bm{\omega},L)]\approx-4.73\times 10^{-6}p_{\text{dir}}+8.35\times 10^{-6}$, respectively. Interestingly, while the both slopes are negatively sloped, only the expected value $\mathbb{E}[J(\bm{\omega},L)]$ displays a sufficiently strong trend, even though significant fluctuations are present. Thus, increasing the directedness of a network appears to have a significant impact on the expected value of the SAF, but relatively little effect on the variance.

4.3 Degree-Frequency Correlations Promote Synchronization

The last property we investigate is the presence of correlations between local networks structure and natural frequencies, i.e., structure-dynamics assortivity [18]. Previous work has shown that positive correlations between oscillators’ nodal degrees and natural frequencies can give rise to explosive synchronization and other novel dynamical behaviors [9, 35, 30]. Moreover, positive correlations between oscillators’ nodal degree and the absolute value of natural frequencies appear in oscillator networks engineered for optimal synchronization [5, 20, 23, 36, 37]. However, especially in this latter case, the observation of correlations in optimized oscillator networks does not imply that such correlations generically promote synchronization properties. Here, we use the framework of Section 3 to investigate this specific relationship.

From the viewpoint of uncertain frequencies developed herein, correlations between network structures and specific instances of natural frequencies may arise from properties related to either the means or the variances or both of these latent variables. We focus on two possible kinds of correlations: (i) correlations between the oscillators’ nodal degrees $k_{i}$ and natural frequency magnitudes $|\mu_{i}|$ and (ii) correlations between the oscillators’ nodal degrees $k_{i}$ and the natural frequency variances $\sigma_{i}^{2}$. For a given oscillator network structure such correlations may be measured using the Pearson correlation coefficient

[TABLE]

for generic nodal properties $x_{i}$ and $y_{i}$ for $i=1,\dots,N$, so that degree-frequency mean and degree-frequency variance correlations may be measured using $\rho_{k,|\mu|}$ and $\rho_{k,\sigma^{2}}$, respectively.

For simplicity we consider two different cases. To study degree-frequency mean correlations we let means $\mu_{i}$ be drawn randomly and uniformly from the interval $[-5,5]$ and natural frequencies be drawn independently (so that covariances are zero) with equal variance $\mathrm{Var}(\omega)$, in which case $\mathbb{E}[J(\bm{\omega,L})]$ and $\mathrm{Var}[J(\bm{\omega},L)]$ are described by Corollary 3.7 and Eq. 17, respectively. To study degree-frequency variance correlations we let means $\mu_{i}=\mu=0$ be identically zero and natural frequencies be drawn independently (so that covariances are zero) with variances $\sigma_{i}^{2}$ drawn randomly and uniformly from the interval $[1,10]$, in which case $\mathbb{E}[J(\bm{\omega,L})]$ and $\mathrm{Var}[J(\bm{\omega},L)]$ are described by Corollary 3.5 and Eq. 16, respectively. We investigate the effect that degree-frequency mean and degree-frequency variance correlations have in collection of ER networks of size $N=500$ with mean degree $\langle 10\rangle$. For each network, after means $\mu_{i}$ or variances $\sigma_{i}^{2}$ are initially drawn, we rearrange them on the network to attain a correlation coefficient $\rho$ randomly drawn from the interval $[-1,1]$.

In Fig. 11, we present the result from degree-frequency mean correlations, plotting the expected value $\mathbb{E}[J(\bm{\omega},L)]$ and variance $\mathrm{Var}[J(\bm{\omega},L)]$ of the SAF versus the correlation coefficient $\rho_{k,|\mu|}$ in panels (a) and (b), respectively. We note a strong relationship between $\mathbb{E}[J(\bm{\omega},L)]$ and $\rho_{k,|\mu|}$, captured roughly by the line $\mathbb{E}[J(\bm{\omega},L)]\approx-1.07\times 10^{-1}\rho_{k,|\mu|}+2.03\times 10^{-1}$, indicating that positive correlations between degrees $k$ and natural frequency mean magnitudes $\mu$ significantly promotes synchronization. On the other hand, a weaker relationship exists between $\mathrm{Var}[J(\bm{\omega},L)]$ and $\rho_{k,|\mu|}$, captured roughly by the line $\mathrm{Var}[J(\bm{\omega},L)]\approx-1.36\times 10^{-4}\rho_{k,|\mu|}+1.55\times 10^{-4}$, indicating that positive correlations between degrees $k$ and natural frequency mean magnitudes $\mu$ decrease the variance of the SAF, but not significantly.

Overall, we find a similar results for the case of degree-frequency variance correlations, presented in Fig. 11, where we plot the expected value $\mathbb{E}[J(\bm{\omega},L)]$ and variance $\mathrm{Var}[J(\bm{\omega},L)]$ of the SAF versus the correlation coefficient $\rho_{k,\sigma^{2}}$ in panels (a) and (b), respectively. Again, we note a strong relationship between $\mathbb{E}[J(\bm{\omega},L)]$ and $\rho_{k,\sigma^{2}}$, captured roughly by the line $\mathbb{E}[J(\bm{\omega},L)]\approx-4.07\times 10^{-2}\rho_{k,\sigma^{2}}+1.17\times 10^{-1}$, indicating that positive correlations between degrees $k$ and natural frequency mean variances $\sigma^{2}$ significantly promotes synchronization. On the other hand, a weaker relationship exists between $\mathrm{Var}[J(\bm{\omega},L)]$ and $\rho_{k,\sigma^{2}}$, captured roughly by the line $\mathrm{Var}[J(\bm{\omega},L)]\approx-3.49\times 10^{-4}\rho_{k,\sigma^{2}}+3.30\times 10^{-4}$, indicating that, again, positive correlations between degrees $k$ and natural frequency mean variances $\sigma^{2}$ decrease the variance of the SAF, but not significantly.

5 Discussion

In this paper, we have investigated the behavior of the synchrony alignment function, or SAF, (see Section 2.2 for a definition), which is a framework for evaluating how the interplay of structural network properties and internal dynamics (i.e., natural frequencies) shape the overall synchronization properties of heterogeneous oscillator networks. Extending previous research of the SAF [36, 37, 38, 37, 43, 37] that treats the oscillators’ natural frequencies as deterministic parameters, in Section 3 we considered the oscillators’ natural frequencies as random variables and derived analytical expressions for the expectation and variance of the SAF. Theorems 3.3 and 3.13 and their corollaries characterize the expectation and variance of the SAF for when the frequencies are uncertain—they are drawn from a distribution with given mean and variance. Knowledge of the precise distribution is not required.

An important insight from Theorem 3.3 is that the expected SAF yields two terms (one depends on the mean frequencies $\bm{\mu}$ and the other depends on the covariance matrix $\Sigma$), and one can consider optimizing these two terms independently. The first term is the SAF, $J(\bm{\mu},L)$, with the expected frequencies as the first argument, rather than the particular frequencies $\bm{\omega}$. This provides theoretical support for why one can optimize the SAF using approximated frequencies, even when the exact frequencies are unknown [see also figure 6.2(c) in [43] for a related numerical experiment]. One can optimize $J(\bm{\mu},L)$ under a variety of system constraints on $L$ and/or $\bm{\mu}$ using the same techniques that were explored in previous research on the SAF [36, 37, 43]. For example, one can consider fixing the network structure and tuning the oscillators’ frequencies, fixing the oscillators’ frequencies and tuning the network structure design, or tuning both the network structure and the oscillators’ frequencies. Optimizing the second term in Theorem 3.3, $\frac{1}{N}\mathrm{Tr}[U(S^{\dagger})^{2}U^{T}\Sigma]$ is less straightforward and one’s strategy may differ depending on the assumptions about $\Sigma$.

Rather then explore the application of these results for optimization, however, we have applied them to gain fundamental insight into which network structure and dynamical properties promote synchronization. In particular, in Section 4 we provided theoretical and experimental for how link weight delocalization, link directedness, and degree-frequency correlations all enhance the synchronization properties of heterogeneous oscillator systems. For the case of IID-distributed random frequencies, we showed the expectation and variance of the SAF are proportional to the second and fourth moments of the spectral density (singular values in the case of directed networks and eigenvalues in the case of undirected networks). Intuitively, these spectral moments can be interpreted as measures for spectral heterogeneity, and our framework predicts synchronization is promoted by spectral concentration, i.e., a lack of spectral heterogeneity. This finding establishes an interesting connection to the study of synchronization for identical oscillators [3, 19, 21, 42], for which the eigenratio — one measure for spectral heterogeneity — is a widely accepted measure for the synchronizability of a system [3]. Our findings for heterogeneous-but-uncertain frequencies are therefore in agreement with, but notably different from this theory. We expect many applications can be benefited by treating the oscillator frequencies as uncertain — that is rather than identical or deterministic — and the results presented herein provide an important step in this direction.

Appendix A Proof of Theorem Theorem 3.1

Proof A.1**.**

This result is well known and is given as Theorem 5.2a in [25]. To provide the reader insight, we provide a simple proof here:

[TABLE]

Appendix B Proof of Theorem Theorem 4.2

We begin by considering, not the sequence of matrices $L^{(m)}$, but $\widehat{L}^{(m)}=L^{(m)}/w$. Note that the new matrices $\widehat{L}^{(m)}$ have off-diagonal entries that are either [math] or $-1$, and therefore can be treated with Theorem 4.1. However, the new node degrees associated with the matrix $\widehat{L}^{(m)}$ is $k/\ell$, so the spectral densities limit to the distribution

[TABLE]

Now, since each $L^{(m)}=\ell\widehat{L}^{(m)}$ we have that $\widehat{\lambda}_{n}^{(m)}=\lambda_{n}^{(m)}/\ell$. Thus, we have that as $m\to\infty$ the spectral density functions $\rho_{m}(\lambda)\to\rho(\lambda)$, where

[TABLE]

which completes the proof.