Chimera states in neural networks and power systems
Partial, frustrated synchronization and chimera-like states are expected to
occur in Kuramoto-like models if the spectral dimension of the underlying graph
is low: $d_s < 4$. We provide numerical evidence that this really happens in
case of the high-voltage power grid of Europe ($d_s < 2$), a large human
connectome (KKI113) and in case of the largest, exactly known brain network
corresponding to the fruit-fly (FF) connectome ($d_s < 4$), even though their
graph dimensions are much higher, i.e.: $d^{EU}_g\simeq 2.6(1)$ and
$d^{FF}_g\simeq 5.4(1)$, $d^{\mathrm{KKI113}}_g\simeq 3.4(1)$. We provide local
synchronization results of the first- and second-order (Shinomoto) Kuramoto
models by numerical solutions on the FF and the European power-grid graphs,
respectively, and show the emergence of \red{chimera-like} patterns on the
graph community level as well as by the local order parameters.
Shengfeng Deng, Géza Ódor
# Chimera states in neural networks and power systems

###### Abstract

Partial, frustrated synchronization and chimera states are expected to occur in Kuramoto-like models if the spectral dimension of the underlying graph is low: \(d_{s}<4\). We provide numerical evidence that this really happens in case of the high-voltage power grid of Europe (\(d_{s}<2\)) and in case of the largest, exactly known brain network corresponding to the fruit-fly (FF) connectome (\(d_{s}<4\)), even though their graph dimensions are much higher, i.e.: \(d_{g}^{EU}\simeq 2.6(1)\) and \(d_{g}^{FF}\simeq 5.4(1)\), \(d_{g}^{KKII113}\simeq 3.4(1)\). We provide local synchronization results of the first- and second-order (Shinomoto) Kuramoto models by numerical solutions on the the FF and the European power-grid graphs, respectively, and show the emergence of chimera-like patterns on the graph community level as well as by the local order parameters.

**We show that Kuramoto oscillator models on large neural connectome graph of the fruit-fly, a human brain, as well as on the power-grid of Europe produce chimera states. This is in agreement with the low spectral dimensions that we calculated by the eigenvalue spectra of the Laplacian of these networks. We compare these results with the topological dimension measurements and previous simulations, strengthening that frustrated synchronization should occur, which can generate slow relaxations, obtained in previous studies within the neighborhood of the synchronization transition point.**

+
Footnote †: preprint: APS/123-QED

## I Introduction

Synchronization phenomena are very widespread in nature and the understanding of their behavior is in the focus of interest. In neural systems, like the brain oscillatory behavior of building elements has been measured by different techniques, while in case of power grids, the alternating currents can also be described by coupled oscillators. Both systems are expected to operate close to the synchronization transition point. In case of the normal brain, self-tuning to the critical point is hypothesized [1] and confirmed by experiments [2] and theoretical considerations [3]. The advantage of criticality is the optimal computational performance, sensitivity as well as dynamically generated long-range memory and interactions [4]. In case of power grids the competition of supply and demands tune the system close to the synchronization transition point [5].

Synchronization models described by the first Kuramoto equation [6] have recently been investigated on complex networks and partial synchronization was found if the spectral dimension is below 4 even if generalizations of the Euclidean dimension, the graph and the Hausdorff dimension are high or diverge [7]. Partial synchronization is more probable in strongly connected modules or communities, which also happens both in biological and technical structures. Modular and most often hierarchical organization is known in general brain networks [8], among others in case of the fruit-fly (FF) connectome [9; 10], as well as in power grids [11]. Thus synchronization occurs in the strongly coupled modules first, while in the loosely coupled parts, nodes may remain desynchronized for the same conditions, which was called frustrated synchronization [12], reminiscent of the semi-critical Griffiths Phases (GP) of condensed matters [13]. Besides, in these phenomena, fluctuations of the global order parameters diverge in an extended control parameter space. Recently this was shown in case of Kuramoto models on brain connectomes [14; 15; 10; 12; 16] as well as in case of power-grids [17; 18; 11].

One can also relate such structures, emerging in these heterogeneous systems, to chimera states, in which subsets of an ensemble of identical, interacting oscillators exhibit distinct dynamical states, such as one group of synchronized oscillators and one group of desynchronized oscillators [19]. Firstly chimeras were defined in systems of identical oscillators [19; 20]. In such a case, a non-zero phase lag term is essential for partial synchronization to occur. Realistic models, however, require oscillators to be heterogeneous and chimeras have been detected on complex [21; 22], brain-like networks [23]. The purpose of the present study is to provide numerical evidences of such structures by solving Kuramoto equations in seemingly different areas of complex systems, on the largest available brain connectome and on the European high voltage power-grid network. We show they are characterized by a low spectral dimension: \(d_{s}<4\). In these models quenched heterogeneity is present structurally, by the topology of the graphs as well as by the different self-frequencies of nodes.

## II Methods

In this section, we detail the models and the methods we applied to describe synchronization on different networks.

### The first-order Kuramoto model

Several oscillator models have been used in biology, the simplest possible one is the Hopf model [24], which has been used frequently in neuroscience, as it can describe a critical point with scale-free avalanches, with sharpened frequency response and enhanced input sensitivity. The local dynamics of each brain area (node) is described by the normal form of
a supercritical Hopf bifurcation, also called a Landau-Stuart oscillator, which is the canonical model for studying the transition from noisy to oscillatory dynamics.

Another complex model, describing more non-linearity [25] is the Kuramoto model [26; 6], with phases \(\theta_{i}(t)\), located at the \(N\) nodes of a network, according to the dynamical equation

\[\dot{\theta}_{i}(t)=\omega_{i}^{0}+K\sum_{j}W_{ij}\sin[\theta_{j}(t)-\theta_{i} (t)]\,. \tag{1}\]

The global coupling \(K\) is the control parameter of this model, by which we can tune the system between asynchronous and synchronous states. The summation is performed over the nearest neighboring nodes, with connections described by the weighted/unweighted adjacency matrix \(W_{ij}\) and \(\omega_{i}^{0}\) denotes the intrinsic frequency of the \(i\)-th oscillator. For simplicity, we used the Gaussian distribution with zero mean and unit variance for the self-frequency distribution \(g(\omega_{i}^{0})\) with respect to a rotating frame [14]. Using this model the resting state critical behavior on large human connectomes [14; 15] was compared with that of the FF [10] on the global order parameter level and the topology dependence has been pointed out, which suggested extended fluctuation region and GP-like behavior in case of the human connectomes in contrast with the FF network.

Very recently we have also investigated an extension of Eq. (1) to the Shinomoto-Kuramoto (SK) model, with periodically driven forces [27] to describe task phase of the brain models [16]

\[\dot{\theta}_{j}(t) = \omega_{j}^{0}+K\sum_{k}W_{jk}\sin[\theta_{k}(t)-\theta_{j}(t)]\] \[+ F\sin(\theta_{j}(t))+\epsilon\eta_{j}(t)\.\]

Here \(\epsilon\) describes an excitation, with a zero centered, Gaussian random annealed noise \(\eta_{j}(t)\) and a site-dependent periodic force term, proportional to a coupling \(F\), was also added. But in fact, a small \(\eta\) proved to be irrelevant, for the synchronization transition, caused by \(F\), in the presence of the chaotic noise. One of the main conclusions of Ref. [16] was that community dependent values of the Hurst exponent \(H\) and the \(\beta\) exponent, measuring self-similarity of time series, varied more with \(F>0\), than in the resting state of the brain, corresponding to \(F=0\). Now we shall test this community dependence of \(R\) and \(\Omega\) in the steady state.

### The second-order Kuramoto model

The time evolution of power-grid synchronization is described by the swing equations [28], set up for mechanical elements (e.g. rotors in generators and motors) with inertia. It is formally equivalent to the second-order Kuramoto equation [29], for a network of \(N\) oscillators with phases \(\theta_{i}(t)\):

\[\dot{\theta}_{i}(t) = \omega_{i}(t) \tag{3}\] \[\dot{\omega}_{i}(t) = \omega_{i}^{0}-\alpha\dot{\theta}_{i}(t)+K\sum_{j=1}^{N}A_{ij}\sin [\theta_{j}(t)-\theta_{i}(t)]\.\]

Here \(\alpha\) is the damping parameter, which describes the power dissipation, or an instantaneous feedback [17], \(K\) is the global coupling, related to the maximum transmitted power between nodes; and \(A_{ij}\), which is the adjacency matrix of the network, contains admittance elements. The quenched external drive, denoted by \(\omega_{i}^{0}\), which is proportional to the self-frequency of the \(i\)-th oscillator and carries a dimension of inverse squared time \([1/s^{2}]\), describes the power in/out of a given node when Eq. (3) is considered to be the swing equation of a coupled AC circuit, but here, similar to the first-order Kuramoto model, we have chosen it zero centered Gaussian random variable as rescaling invariance of the equation allows to transform it out within a rotating frame. For simplicity, one can assume that \(\omega_{i}(0)\) is drawn from the same distribution as \(\omega_{i}^{0}\) and numerically set \(\omega_{i}(0)=\omega_{i}^{0}\), amounting to taking \([s]\)=1. In our present study the following parameter settings were used: the dissipation factor \(\alpha\), is chosen to be equal to 0.4 to meet expectations for power grids, with the \([1/s]\) inverse time physical dimension assumption.

To characterize the phase transition properties the phase order parameter \(R(t)\) has been studied for both the first- and second-order Kuramoto models. To ensure the relaxation to the steady states, we measured the Kuramoto phase order parameter

\[z(t_{k})=r(t_{k})\exp[i\theta(t_{k})]=1/N\sum_{j}\exp\left[i\theta_{j}(t_{k}) \right]\, \tag{4}\]

where \(0\leq r(t_{k})\leq 1\) gauges the overall coherence and \(\theta(t_{k})\) is the average phase at discrete sampling times \(t_{k}\), which was chosen to follow an exponential growth: \(t_{k}=1+1.08^{k}\) to spare memory space. The calculation of derivatives was done adaptively at small time steps via the Bulirsch-Stoer stepper [30]. The sets of equations (1), (2) and (3) were solved numerically for \(10^{3}-10^{4}\) independent initial conditions in [18], initialized by different \(\omega_{i}^{0}\)-s and different \(\theta_{i}(0)\)-s if disordered initial phases were invoked. Then sample averages for the phases and the frequencies give rise to the Kuramoto order parameter

\[R(t_{k})=\langle r(t_{k})\rangle\, \tag{5}\]

and the variance of the frequencies

\[\Omega(t)=\frac{1}{N}\sum_{j=1}^{N}(\overline{\omega}(t)-\omega_{j}^{2}(t)). \tag{6}\]

We don't discuss the time-dependent behavior of the global order parameters as this has been investigated in detail in [31; 18]. In the steady state, which we determined by visual inspection of \(R(t)\) and \(\Omega(t)\), we measured their half values and the standard deviations \(\sigma(R(t))\) and \(\sigma(\Omega(t))\) in order to locate the transition points. In the paper we used the \(\sigma(R)\), \(\sigma(\Omega)\) values, obtained by sample and time averages in the steady state.

### Topological and spectral dimensions

The effective graph (topological) dimension \(d_{g}\) is defined by

\[N(r)\sim r^{d_{g}}, \tag{7}\]
where we counted the number of nodes \(N(r)\) with chemical distance \(r\) or less from randomly selected seeds and calculated averages over many trials [32; 33]. In most cases, \(d_{g}\) as obtained from this cluster growing method can be served as an estimation for the more rigorously defined Hausdorff dimension \(d_{H}\approx d_{g}\)[34]. In Ref. [10] the graph dimension of the FF was estimated to be \(d_{g}^{FF}=5.4(1)\), while in [18] we provided a value for the unweighted European power-grid \(d_{g}^{EU}=2.6(1)\).

For regular Euclidean lattices, it has been shown that a true transition to global synchronization under the thermodynamics limit is only possible for \(d>4\) in both the first- and the second-order Kuramoto models [35; 36], while for \(d\leq 4\), there is only a crossover from the asynchronous phase to a partial synchronous phase characterized by an increasingly broadened variance of \(R\) and a shifting crossover point \(K_{c}^{\prime}\) as the system size increases [37; 36]. The natural question is then if it is the topological dimension \(d_{g}\) defined in Eq. (7), or equivalently the Hausdorff dimension \(d_{H}\), that dictates the synchronization properties for general networks which may assume non-integer dimensions. Refs. [37; 7] suggested that the synchronization properties of a general network should be related to the spectral dimension derived from the eigenvalue spectrum of the graph Laplacian matrix, and even more so for the so-called complex network manifolds studied therein, which are constructed out of finite-dimensional simplicies but are characterized by an infinite Hausdorff dimension due to their small-world properties [37; 38].

Graph spectral properties of complex networks have been shown to be particularly relevant to the structures of networks [39]. Following Refs. [37; 7], we adopt the normalized Laplacian \(L\) with elements

\[L_{ij}=\delta_{ij}-A_{ij}/k_{i} \tag{8}\]

for unweighted networks, where \(k_{i}\) denotes the degree of node \(i\). Similarly, for weighted networks, the elements of the normalized Laplacian are given by

\[L_{ij}=\delta_{ij}-W_{ij}/k_{i}^{\prime}, \tag{9}\]

where \(k_{i}^{\prime}=\sum_{j}W_{ji}\) denotes the weighted in-degree of node \(i\). The normalized Laplacian has real eigenvalues \(0=\lambda_{1}\leq\lambda_{2}\leq\cdots\leq\lambda_{N}\), the density of which scales as [37; 40]

\[\rho(\lambda)\simeq\lambda^{d_{g}/2-1} \tag{10}\]

for \(\lambda\ll 1\), where \(d_{s}\) is the spectral dimension. The cumulative density is then given by

\[\rho_{c}(\lambda)=\int_{0}^{\lambda}d\lambda^{\prime}\rho(\lambda^{\prime}) \simeq\lambda^{d_{g}/2}. \tag{11}\]

Since Eqs. (10) and (11) hold for small \(\lambda\) values, for the fruit-fly connectome and the European high-voltage power grid network that are going to be studied in Sec. III, with which \(N\gg 1\), we will only extract the densities for the first 200 smallest eigenvalues for ease of eigenvalue computation without loss of generality.

As illustrated in Fig. 1, Euclidean lattices in dimension \(d\) have spectral dimension \(d_{s}=d\). Therefore in this case the spectral dimension is also equal to the Hausdorff dimension of the lattice, \(d_{s}=d_{H}\). However, in general, networks can have non-integer spectral dimension \(d_{s}\) not equal to their Hausdorff dimension. Ref. [37] demonstrated that in lower spectral dimensions \(d_{s}<4\), there is a parameter regime that exhibits frustrated synchronization with spatio-temporal fluctuations even in the stationary state. Then, similar to the emergence of rare regions in Griffiths phases [13], one should expect to observe states with rare regions-usually called "chimera states" - in such frustrated synchronization as well, as we will demonstrate in what follows with the aid of the local order parameter of Kuramoto models.

### Local order parameters of the first- and second-order Kuramoto models

To investigate the heterogeneity further, we measured the local Kuramoto order parameter, defined as the partial sum of phases for the neighbors of node \(i\)

\[r_{i}(t)=\frac{1}{N_{\text{i.neigh}}}\left|\sum_{j}^{|N_{\text{i.neigh}}}A_{ij }e^{i\theta_{j}(t)}\right|. \tag{12}\]

This local Kuramoto measure was firstly suggested by Restrepo _et al._[41; 42] to quantify the local synchronization of nodes, which allows us to visualize regions of synchronized/unsynchronized chimera-like behavior and which will be the main quantity of interest of this paper.

## III Chimera states

### Chimera states in the fruit-fly (FF) and in a human connectome

First, we examine chimera states in the first-order Kuramoto model on the FF connectome. Connectomes are

Figure 1: Spectral dimensions of \(d\)-dimensional regular lattices (\(d=1,2,3,4\)) of different lateral sizes \(l\), measured through Eq. (11). The real dimension of a lattice is asymptotically approached as \(l\rightarrow\infty\).
defined as structural networks of neural connections of the brain [8]. For the fruit fly, we used the hemibrain dataset (v1.0.1) from [43], which has \(N_{FF}=21\,662\) nodes and \(E_{FF}=3\,413\,160\) edges, out of which the largest single connected component contains \(N=21\,615\) and \(E=3\,410\,247\) directed and weighted edges, with weights being the number of connections between a pair of nodes. The number of incoming edges varies between 1 and 2708. The weights are integer numbers, varying between 1 and 4299. The average node degree is \(\langle k\rangle=315.129\) (for the in-degrees it is: 157.6), while the average weighted degree is \(\langle w\rangle=628\). The adjacency matrix, visualized in [10], shows a weak hierarchical modular structure, however, it is not random. For example, the degree distribution is much wider than that of a random graph and exhibits a fat tail. The analysis in [10] found a weight distribution \(p(W_{ij})\) with a heavy tail, and assuming a power-law (PL) form, a decay exponent \(2.9(2)\) could be fitted for the \(W_{ij}>100\) region. The modularity quotient of a network is defined by [44]

\[Q=\frac{1}{N\langle k\rangle}\sum_{ij}\left(A_{ij}-\frac{k_{i}k_{j}}{N\langle k \rangle}\right)\delta(g_{i},g_{j})\,. \tag{13}\]

The maximum of this value corresponds to the optimal community structure characterizes how modular a network is, where \(\delta(g_{i},g_{j})\) is 1 when nodes \(i\) and \(j\) were found to be in the same community \(g\), or 0 otherwise. Community detection algorithms based on modularity optimization get the closest to the actual modular properties of the network. The modularity was calculated using community structures detected by the Louvain method [45], from which we obtained \(Q_{FF}\approx 0.631\)[16]. The effective graph (topological) dimension, obtained by the breadth-first search algorithm is \(d_{g}^{FF}=5.4(5)\).

To compute the spectra dimensions, we extracted the cumulative density distributions of the first 200 eigenvalues of both the Laplacian matrices of the unweighted and weighted FF connectomes and plot them in a log-log scale as shown in Fig. 2. For small enough \(\lambda\) values, the distributions indeed display a scaling regime, permitting estimation of spectral dimensions by Eq. (11), as listed in the plot legend and Table 1. Now, even though \(d_{g}^{FF}>4\), since \(d_{s}<4\), one should expect frustrated synchronization in some parameter regime where chimera states may be observed.

To that end, we solved Eq. (1) and Eq. (2) numerically with respect to different coupling strengths and \(F\) on the weighted FF connectome and observed the respective global phase order parameter \(R(t\rightarrow\infty)\) in the stationary state. The stationary state is typically reached after a few hundred time steps; see, for example, the first panel of Fig. 3. Practically, we followed the dynamics up to \(t=1000\) to ensure stationarity. In Ref. [10], we had estimated the critical coupling \(K_{c}\simeq 1.6\) from the peak of the variance of \(R(t\rightarrow\infty)\) with respect to \(K\). By using this critical coupling, we further calculated the local order parameter Eq. (12) for each node, averaged over 20 independent simulation runs. In Fig. 3, the local order parameters for three representative time steps are displayed by encoding the respective values to a color map. Since the simulations were started from a fully asynchronous state, we see that the system gradually evolves into a more synchronous state at larger times. However, even in the globally stationary state characterized by a constant \(R\) value, the local order parameters show rather inhomogeneous patterns, with some parts of the connectome are more synchronized (greener regions) while some other parts are less synchronized (redder regions), indicating the emergence of chimera states [41; 42]. Note the disparity of the synchronization levels between different regions is quite large in this case, with greener regions almost fully synchronized and red regions fully unsynchronized. What is more, as partially shown by the second and the third panels at \(t=748\) and \(t=1885\), the distribution of the local order parameters can still evolve in the globally stationary state. Simulation results seemed to suggest quite random temporal behavior for the local order parameters (not shown here), but more careful studies for the long-time behavior are still needed to examine if it is periodic with a very long period. These results are thus suggestive of strong spatio-temporal fluctuations in chimera states, as it is typical for frustrated synchronization [37].

To provide more evidence of the Chimera states, we have calculated the order parameters in the steady state at \(K=1.6\) in the nine largest communities, determined by the Louvain method [16]. However, the community dependence is rather weak in case of the Kuramoto model. So we enhanced the local synchronization by adding periodic forces within the framework of the SK model. The transition point shifts in the range \(F_{c}\in[0.05,0.1]\). Even more evident community dependence could be found in the frequency synchronization points estimated by the peaks of the variances of the order parameter \(\Omega\). As one can see in Fig. 6, frequency entrainment occurs in the range \(F_{c}^{\prime}\in[0.025,0.1]\) in different communities.

\begin{table}
\begin{tabular}{l c c} \hline \hline  & \multicolumn{2}{c}{Fruit-fly conn. European grid} \\ \hline \(d_{g}\) & 5.4(1)[10] & 2.6(1)[18] \\ \(d_{s}\) (unwei.) & 3.375 & 1.717 \\ \(d_{s}\) (wei.) & 2.342 & 1.507 \\ \hline \end{tabular}
\end{table}
Table 1: The graph dimensions and spectral dimensions of the FF connectome and the European grid network. Both spectral dimensions for unweighted and weighted networks are given.

Figure 2: The cumulative eigenvalue densities of the unweighted (blue solid square) and the weighted (red dot) FF network, extracted from the first 200 smallest eigenvalues.
That means that for certain forces some communities are in the super-, while others are in the sub-critical states locally, suggesting Chimeras.

As it has already been shown in Refs. [10; 16] the FF graph exhibits a weak modular structure. Much higher level of modularity can be observed in human connectomes, albeit on a coarse grained scale, describing the white matter. The large human connectomes obtained by DTI of MRI [10; 16; 46; 47] has node numbers of order of a million. We have not been able to calculate the local order parameters and the spectral dimensions for such large systems. In [47] we calculated their graph dimensions, which proved to be above 3 of the embedding space, but lower than 4, due to the long fiber tracts, connecting distant regions. We show here, that synchronization of communities of the KKI-113 graph exhibits much more visible differences than that of the FF, suggesting strong Chimeras in case of the Kuramoto model running on them. This small-world graph contains 799 133 nodes, connected via 48 096 500 undirected and weighted edges and exhibit a hierarchical modular structure, because it was constructed from cerebral regions of the Desikan-Killany-Tourville parcellations, which is standard in neuroimaging [48]. The modularity quotient is much higher, than that of the FF: \(Q_{KKI-113}\approx 0.915\) and the topological dimension is just \(d_{g}=3.4(1)\)[16].

As one can see on Fig. 7 in certain communities the synchronization is high at \(K=3\), while others are still practically unsynchronized at this coupling, suggesting Chimera states.

Figure 4: The cumulative eigenvalue densities of the unweighted (red solid square) and the weighted (blue dot) European power-grid network, extracted from the first 200 smallest eigenvalues.

Figure 3: Evolution of \(R(t)\) towards the steady state following a start from a fully unsynchronized state (upper left panel) and local Kuramoto results at five different time steps encoded by the color map of \(r_{i}\) for the weighted FF connectome. Results were averaged over 20 samples. In the color-maps, red corresponds to low local synchronization, and green to high synchronization.
### Chimera states in the European high-voltage power grid network

Unlike neural networks on which the oscillators are massless a power-grid network are massive and should be described by the second-order Kuramoto model. In this subsection, we attempt to show if chimera states can emerge in such systems. Power-grid networks are genuinely hierarchical modular networks if the detailed information for the medium- and low-voltage parts of the grids are also incorporated. Practically, it is almost impossible to infer the entire structure of large power-grid networks, but it is feasible to mimic it by adding medium- and low-voltage parts to the high-voltage (HV) skeleton, according to the empirical hierarchical distribution, as it was done in Ref. [11] We downloaded the European HV power grid from the "SciGRID Dataset"[49] encoding the 2016 status, deduced via processing google street-map. We have not supplemented this graph with lower-voltage parts, but it already contains 12 kV, 20 kV,... links, which belong to

Figure 5: Local Kuramoto results in the stationary state encoded by the color maps of \(1-r_{i}\) for (a) the unweighted and (b) the weighted European high-voltage networks. Results were averaged over 20 samples. Red corresponds to low local synchronization, and green to high synchronization. The width of the gray edges in (b) is proportional to the logarithm of the weights.

Figure 6: Community dependence of \(\sigma(\Omega)\) at \(K=1.6\) in the FF. Peaks, marking the local entrainment point \(F_{c}^{\prime}\), move from 0.025 to 0.1 for communities 8,6,4,3,2,1 and the whole (top to bottom curves). Results were obtained by averaging over 200 samples.

Figure 7: Community dependence of \(R\) for different \(K\)-s showing different phase synchronizations, corresponding to Chimera states in the KKI-113 connectome. The the lowest curve denotes the synchronization of the whole system, which grows the slowest by increasing \(K\). Inset: Fluctuations of the same data, showing peaks at different \(K\)-s. The lowest curve, representing the whole system grows the slowest. Results were averaged over 100 samples.
the middle-voltage category, according to the definition of the 100 kV threshold for HV lines.

This graph contains \(N=13478\) nodes, interconnected via \(E=33844\) links. After symmetrizing it, an average degree \(\langle k\rangle=2.51\) was obtained. In Fig. 1 of Ref. [18] the degree distribution is shown. The tail of the degree distribution for \(k\geq 15\) could be well fitted by a stretched exponential \(8.25\times e^{-0.53(5)k}\) function, which renders this network at the threshold of robust/fragile: \(\gamma=3/2\), as according to the definition in [50], networks with a \(P(k>K)=Ce^{-k/7}\) cumulative degree distribution and \(\gamma<3/2\) are robust, based on a mean-field percolation theory under random node removals. The adjacency matrix, visualized on Fig. 2 of Ref. [18] proves that this is a highly modular graph, characterized by \(Q^{EU}=0.963\). Furthermore, it is a small-world network according to the definition of the small-worldness coefficient [18; 51]. By calculating the graph dimension using the breadth-first search algorithm as shown in the inset of Fig. 3 of Ref. [18], \(d_{g}^{EU}=2.6(1)\) was obtained. Since the coupling between a pair of nodes of a power grid is proportional to the maximal power \(P_{ij}\) transmitted between them and inverse to the imaginary part \(X_{ij}\) of the impedance of the transmission line, weights computed from the normalized values of \(P_{ij}/X_{ij}\) had also been considered to construct the weighted network.

By again extracting the cumulative density distributions of the first 200 eigenvalues of both the Laplacian matrices of the unweighted and weighted European power-grid networks, Fig. 4 shows quite clean power laws. As listed in the plot legend and Table 1, the estimated spectral dimensions are both well below the critical dimension \(d_{c}=4\). Hence, one can again expect to observe chimera states, although instead, the second-order Kuramoto model (3) is going to be inspected in this case.

We again tune the system to the verge of criticality. By solving Eq. (3) to obtain the peak of the variance of \(R(t\rightarrow\infty)\), the critical couplings had been estimated to \(K_{c}=80\) for the unweighted network [18] and \(K_{c}=7000\) for the weighted network. The stationary states are typically reached after a few hundred time steps (see Ref. [18]), but we solved them up to \(t=20000\) to ensure the stationarity of the system. The local order parameters, calculated in the stationary states after averaging over 100 samples, are then obtained with respect to these critical couplings. Fig. 5 shows that inhomogeneous patterns, encoded again in a color map, indeed overwhelm the system. Due to the higher levels of synchronization, there are quite some proportions of oscillators with their local order parameters relatively closer to 1 as compared to the FF connectome case. Hence the color map shows the quantity \(1-r_{i}\) instead. Since the differences in the local parameters in greener regions and redder regions are still quite apparent, we see that as suggested by the low spectral dimension, chimera states indeed can be observed in this case. Note that even though the weighted network is a bit less synchronous globally at \(K=7000\) [\(R(t\rightarrow\infty)\simeq 0.47\)] than the unweighted network at \(K=80\) [\(R(t\rightarrow\infty)\simeq 0.48\)], the weighted network still seems to be more synchronized locally in many parts of the network. This emphasizes the importance of incorporating edge weights to take into account more realistic couplings between the nodes.

Comparing the local order parameter patterns in Fig. 3 and Fig. 5, it is also interesting to note that less synchronous regions are typically also less clustered as compared to regions with higher levels of synchronization. This is in some sense in reminiscence of the analysis in Ref. [12], in which it had been shown that chimera states can also be characterized by the order parameters of different moduli.

To provide more evidence of the Chimera states, we have also calculated the steady-state \(R\) in the twelve largest communities, determined via the Louvain method with a modularity score close to the maximum \(Q\approx 0.795\), in the same way as in case of the FF [52]. As one can see in Fig.8, synchronization occurs at different couplings in different communities, such that for small \(K\)-s the small communities are fully ordered, while the larger ones are still desynchronized. This is related to the size dependence of \(K_{c}\) in case of crossover, however here the the communities are not independent. Note, that the fully ordered communities have less than 100 nodes.

## IV Summary

In this paper, we have demonstrated that chimera states can occur in Kuramoto-type models on large networks if the spectral dimension is low, i.e. \(d_{s}<4\), even if the graph dimension is not necessarily like that. That happened in case of the graph of the FF connectome, which exhibits \(d_{g}=5.4(1)\). This is in agreement with the hypothesis, advanced for the first-order Kuramoto model in [7]. But as modularity is weak for FF, so do

Figure 8: Community dependence of \(R\) for different \(K\)-s showing different phase synchronizations, corresponding to Chimera states. The thick black curve denotes the synchronization of the whole system, which grows the slowest by increasing \(K\). Inset: Fluctuations of the same data, showing different synchronization points. The thick black curve, representing the whole system has the rightmost peak. Results were obtained by averaging over 100 samples.
the Chimeras. We can show them by a community-level analysis with an applied periodic external field. In contrast, for a large human connectome possessing high modularity, we show strong community dependence of the local synchronization.

Power grids can be described by the second-order Kuramoto model, which possesses inertia. We found that the European HV power grid has a graph dimension \(d_{g}=2.6(1)\), but the spectral dimensions seem to be below \(d_{s}=2\). Still, the occurrence of chimera-like patterns can be observed via order parameters and confirmed by a community-level synchronization study. We demonstrated the level of local synchronization by showing the local Kuramoto order parameter, but similar results have been found by calculating the local frequency spreads.

###### Acknowledgements.

 We thank Kristof Benedek and Balint Hartmann for providing weight calculations of the European network, Istvan Papp for exploring the communities, Jeffrey Kelling for developing the GPU solver code, and Robert Juhasz for the helpful discussions. This research was funded by ELKH grant SA-44/2021, and the Hungarian National Research, Development, and Innovation Office NKFIH grant K128989. Most of the numerical work was done on KIFU supercomputers of Hungary.

## Data Availability Statement

Data are available on request from the corresponding author.
