Long Wavelength Coherency in Well Connected Electric Power Networks

We investigate coherent oscillations in large scale transmission power grids, where large groups of generators respond in unison to a distant disturbance. Such long wavelength coherent phenomena are known as inter-area oscillations. Their existence in networks of weakly connected areas is well captured by singular perturbation theory. However, they are also observed in strongly connected networks without time-scale separation, where applying singular perturbation theory is not justified. We show that the occurrence of these oscillations is actually generic. Applying matrix perturbation theory, we show that, because these modes have the lowest oscillation frequencies of the system, they are only moderately sensitive to increased network connectivity between well chosen, initially weakly connected areas, and that their general structure remains the same, regardless of the strength of the inter-area coupling. This is qualitatively understood by bringing together the standard singular perturbation theory and Courant's nodal domain theorem.


I. INTRODUCTION
Synchronous generators in interconnected AC electric power systems exhibit electro-mechanical oscillations. Of particular interest are large-scale cooperative phenomena termed inter-area oscillations which are coherent, sub-Hz frequency oscillations between geographically separated large groups of generators [1]. Such oscillations may become unstable and lead to large scale blackouts [2], therefore safe system operation requires that they are appropriately damped, which becomes harder as the energy transition unfolds. As a matter of fact, power system stabilizers installed on conventional synchronous generators have so far been the main source of damping against inter-area oscillations [1] and substituting new renewable sources of energy for conventional synchronous machines reduces the availability of these resources. There is a vast literature on damping of inter-area oscillations in power systems with large penetrations of new renewable generation, see e.g. [3].
To achieve optimal damping of inter-area oscillations it is important to first identify the geographical areas that carry them, i.e., where synchronous generators display the same frequency response following a fault or other excitations. This identification commonly proceeds through highlighting weak links [4], spectral analysis [5], and data-based approaches identifying generators with similar frequency responses, either in simulations of the linearized dynamics [6] or using wide area measurement data [7]. Once coherent areas are identified, inter-area oscillations are studied using aggregated models constructed from singular perturbation theory [5], [8], [9]. These methods presuppose the existence of at least one small parameter µ measuring the ratio between the inter-area and the intra-area connection strength or the associated time scales. When µ is not small, the theory loses its validity, yet inter-area oscillations are observed even in networks with large µ. This is illustrated in Fig. 1 which shows coherent, low-frequency inter-area oscillations obtained numerically in the PanTaGruEl model of the synchronous grid of continental Europe [10], [11]. It is seen that the Iberian Peninsula responds coherently to a noisy power injection in Greece, and reciprocally, even though µ is close to one hundred. All 981 nodes in the Iberian Peninsula respond coherently -with the same frequency and phase -to a single-node, noisy perturbation in the Balkans, as do all 368 Balkan nodes to a similar perturbation in the Iberian Peninsula. The system of Fig. 1 operates well outside the regime of validity of the standard singular perturbation theory [5], [8], [9]. Another viewpoint predicts the existence of inter-area oscillations, regardless of the network connectivity. The dynamics of voltage angles in power systems is commonly modeled by the swing equations, which are a set of coupled ordinary differential equations [12]. The coupling is determined by a graph Laplacian matrix, whose eigenvectors naturally define coherent areas. As a matter of fact, Courant's nodal domain theorem states that the k th eigenfunction of an elliptic operator acting on a bounded domain Ω defines no more than k nodal subdomains of Ω, where the eigenfunction does not change sign [13]. To make a long story short, the eigenmodes with lowest eigenvalues of an elliptic operator such as, e.g., a Laplacian operator, define few large areasnodal subdomains -on which the sign of their components does not change. When a single or very few eigenmodes are excited, the resulting oscillations appear coherent inside the corresponding nodal domains. The theorem has recently been extended from continuous elliptic operators to graphs represented by, e.g., a discrete Laplacian matrix [14]. In the case of power systems, that Laplacian matrix represents a quasi-planar graph and nodal domains are two-dimensional areas. Courant's nodal domain theorem states that these areas are larger for slower eigenvectors of the graph Laplacian -those with lower frequencies. Consequently, these modes do not resolve inhomogeneities in the inertia and damping parameters, therefore the structure of the system's true interarea modes is directly determined by the slowest eigenvectors of the Laplacian. Low frequency coherent oscillations over large areas thus naturally emerge from a modal decomposition of the graph Laplacian. Our purpose in this manuscript is to connect the standard singular perturbation theory approach to inter-area oscillations to this modal point of view, valid regardless of inter-area coupling. We will use matrix perturbation theory [15] to describe low frequency oscillations, extending our earlier work [16]. We will show how an appropriate choice of areas gives an excellent approximation for the modes responsible for interarea oscillations, which is only very weakly sensitive to the inter-area coupling. This fills an important gap in the theory of inter-area oscillations as it explains their persistence in strongly connected power networks.
The manuscript is organized as follows. Section II introduces the network model and motivates the focus on the network Laplacian to describe the inter-area oscillations. Matrix perturbation theory is discussed in Section III. After a summary of the general results of perturbation theory we define our model for describing inter-area oscillations. In the rest of the section the validity of perturbation theory is discussed in terms of series convergence and the change of eigenvectors due to avoided crossings. Applications of the theory to a synthetic two-area network, the IEEE RTS 96 Test System, and the PanTaGruEl model of the synchronous very high voltage grid of continental Europe are shown in Section IV. Results are discussed in Section V.

A. Swing Equations
We use the structure-preserving model of Ref. [17] and consider the voltage angle dynamics of a high voltage power grid with N nodes. The dynamics of the voltage angle θ i on generator nodes is determined by the swing equations [12]. In the case of high voltage transmission grids, a standard approximation is the lossless line approximation, which neglects Ohmic losses. The swing equations then read with the inertia m i and damping d i parameters of the generator and their active power P i > 0. Loads are assumed frequencydependent, and the power they draw is given by a constant P i < 0 and a frequency-dependent term, d iθi . Load voltage angles then obey [17] where the frequency dependence of the loads is determined by the parameter [17] with P 0 i the power consumption at the nominal frequency ω 0 . The parameter α has been evaluated experimentally, α ∈ [0. 8,2] [18], [19]. In Eqs. (1) and (2), B ij denotes the product of the voltage magnitudes at nodes i and j with the line susceptance. In the lossless line approximation, line conductances are neglected.
We investigate the generation of inter-area oscillations through small signal stability analysis. Accordingly, we linearize Eqs. (1) and (2) about the operational synchronous state with θ i = θ where we grouped the voltage angle deviations into a vector δθ, and introduced the diagonal inertia and damping matrices, as well as the network Laplacian matrix L, Eq. (4) is often written as a linear system of ordinary differential equations with stability matrix A, where angles, frequencies and power injections are grouped into x and Π as which defines the A-matrix Above, 0 denotes either the matrix or the vector with all components equal to zero, and I is the identity matrix, subindices g and l denote generators and loads respectively. Inter-area oscillations correspond to the slowest eigenmodes of A. From Eqs. (4) and (6), they are determined by the Laplacian, the inertia, and the damping matrices.

B. Eigenvectors and Eigenvalues
Eigenvectors and -values of A are easily obtained in the case of homogeneous damping d i = d and inertia m i = m. For each mode, the oscillation frequency reads where λ α is an eigenvalue of L and γ = d/m [20]. Furthermore, both the angle and frequency components of each eigenmode of A are given by those of an eigenmode of L. Thus, the oscillation frequency and the mode structure are related to the eigenvalues and -vectors of L, when damping and inertia are homogeneous. The homogeneity assumption is not present in real electric power grids, not even after a Kron reduction of the inertialess load nodes. A question of central interest is therefore how much do inertia and damping inhomogeneities affect the structure of the eigenmodes of A. Our matrix perturbation approach to be presented below shows that long wavelength, slow modes are only weakly affected by such inhomogeneities. We illustrate this numerically on PanTaGruEl in Fig. 2. We write damping and inertia as where the bar indicates the average over the nodes in the network, δx i = x i −x with the real value x i = m i or d i in the system, and χ ∈ [0, 1] is a continuous parameter tuning the system from the homogeneous configuration (χ = 0) to its real configuration (χ = 1). The top panel of Fig. 2 shows that the oscillation frequencies of the slowest modes of PanTaGruEl are only weakly sensitive to inertia and damping inhomogeneities. The overlap |ψ α (χ)ψ α (0)| of the lowest five right eigenvectors at χ = 0 and χ ∈ [0, 1] is next shown in the bottom panel of Fig. 2. The data illustrate nicely that inhomogeneities have only a minor effect on the oscillation frequency and the spatial structure of the lowestlying eigenvectors -those of interest in this manuscript. This can be explained by the long-wavelength nature of these modes, which effectively averages damping and inertia over areas that are so large that the modes do not resolve their inhomogeneities. Ref. [21] found that the change of ω α with χ is to leading order proportional to i u 2 αi δγ i , where u α is the eigenvector of the Laplacian corresponding to λ α and δγ i = (d i /m i − γ). For the lowest modes which are approximately constant over large areas, u 2 αi can be factored out and i∈area δγ i ≈ 0, assuming that the distribution of δγ i is independent of the geographical location. A similar argument holds for the eigenvectors. The corrections to the right eigenvector ψ α of the stability matrix are given by a sum over all the other eigenvectors ψ β , β = α, each weighted by a factor proportional to i u βi u αi δγ i . Assuming again that the lowest modes are mostly constant over large areas, the contributions from the highly fluctuating higher modes average out.
Higher up in the spectrum, modes have a shorter wavelength and hence better geographical resolution, and these statements obviously break down. This is of no consequence for the validity of our approach to slow, large-wavelength modes.

III. MATRIX PERTURBATION THEORY
Consider a matrix, which we are able to diagonalize exactly. Suppose we perturb that matrix as L 0 → L(ε) = L 0 + εL 1 , with a perturbation matrix L 1 that does not commute with L 0 . Matrix perturbation theory [22] is a method for describing the change of eigenvalues and -vectors of L(ε) as a series expansion in the dimensionless parameter ε. It is a standard method of theoretical physics [23] that has recently been exported to electric power and other network systems. It has been used to investigate the change in oscillation frequencies under small changes of certain slow modes in [24]. It has been used in Ref. [25] to construct a control scheme for the output of generators and enhance power grid stability. The optimal placement of inertia and damping has further been investigated using perturbation theory in Ref. [21]. In a more general context, Ref. [26] uses matrix perturbation theory to investigate a network of networks.
Taken as a whole, perturbation theory is valid as long as ε 1. However, we argue below that, when applied to a restricted range of low-frequency modes -such as few of the slowest modes represented in color in Fig. 2 -the validity range generally becomes significantly larger and may even include the ε → 1 limit. Before we apply it to our problem, we first give a brief general description of non-degenerate and degenerate perturbation theory in the next paragraph..
A. General Framework 1) Non-degenerate perturbation theory: Take a real symmetric matrix L 0 with known eigenvalues λ α . This matrix is subjected to a perturbation εL 1 , where L 1 is a symmetric matrix which does not commute with L 0 and ε is a dimensionless scalar parameter. Perturbation theory expands the eigenvalues λ α and -vectors u α of L(ε) = L 0 + εL 1 in a power series in ε, α ) the nth order correction to the eigenvalue (-vector). The first and second order corrections of the eigenvalues read while the first order correction of the eigenvectors reads Higher order corrections can be obtained recursively from the eigenvalue problem using the series expansion (11). For more details, the reader is referred to Refs. [15], [22], [23].
2) Degenerate perturbation theory: Eqs. (12) and (13) are only valid as long as the α th unperturbed eigenvalue λ (0) α has multiplicity one. We call this the nondegenerate case. Special care needs to be taken when considering corrections to eigenvalues with multiplicity larger than one. In this case the corresponding eigenvectors are not unique and span the degenerate subspace D. This degenerate subspace has to be considered separately from the rest of the vector space. The eigenbasis spanning D is a priori not uniquely defined, since any normalized linear combination of the degenerate eigenvectors is also an eigenvector. However there is one and only one linear combination for which the change in eigenvectors is smooth as the perturbation is turned on. Degenerate perturbation theory dictates to choose that linear combination as a starting point. It diagonalizes L 1 within D and is defined by the conditions The first condition in Eq. (15) readily gives the first-order correction to the degenerate set of eigenvalues. Higher-order corrections to eigenvalues and eigenvectors of D are given by Eqs. (12) and (13), with the substitution β = α → β / ∈ D.

B. Specific Set-up
We consider a network partitioned into p, initially disconnected areas labeled a, each containing n a nodes and represented by a n a × n a Laplacian matrix L a . The unper- Because each area is connected, to each block a corresponds a single eigenvalue λ (0) a = 0, associated to an eigenvector being constant in area a and zero everywhere elsẽ where 1 is the vector of ones, and the tilde means that these eigenvectors do not satisfy (15) yet. From now on, we refer to these eigenvectors as zero-modes. They correspond to each area oscillating coherently on its own, and we are interested in finding out how they change when the inter-area coupling is turned on. If they do not change significantly, the area will engage in coherent oscillations, and we will see that this is the case for well chosen areas. We therefore focus on the set D spanned by the zero-modes from here on. The perturbation εL 1 contains the lines that connect the different areas and ε tunes the network from having unconnected areas at ε = 0 to recovering the real, fully connected network at ε = 1. To obtain the linear combination of eigenvectors that satisfy (15) we project L 1 onto D and linearize it. The projected matrix L proj is given by where B ab = i∈a;j∈b B ij is the sum of all connections between area a and b. L proj has a zero eigenvalue with eigenvector The corresponding linear combination then gives a v 0aũa = 1 N / √ N , i.e., the global zero-mode of the full network. The p − 1 other eigenvectors of L proj define p − 1 other linear combinations of zero-modes constituting the unperturbed basis in which the perturbation expansions are constructed. We call these linear combinations "hybridized zero-modes". Our theory to be presented below focuses on them and on how they evolve as the inter-area connections increase.
The hybridized zero-modes acquire nonzero first-order eigenvalues which are linear in ε, with a slope determined by (15). Second-order corrections to their eigenvalues emerge due to the interaction with non-zero-modes triggered by εL 1 . These corrections read Because λ (0) α = 0 for the hybridized zero-modes, the secondorder corrections are in particular negative, reflecting the eigenvalue repulsion [27] between the hybridized and the non-zero-modes. Simultaneously, Gershgorin's circle theorem guarantees that eigenvalues of L 0 + εL 1 are nonnegative [22]. These two effects result in the behavior of the hybridized eigenvalues shown in Fig. 3. The eigenvalues of the hybridized modes are shown in color. We see that after a short linear rise captured by first-order perturbation theory, they all quickly bend downward to reach what looks like a horizontal asymptotic. Furthermore, only the upper one (dark blue curve) gets close to the non-zero modes (gray curves) as ε increases. In the following paragraphs we analyze this behavior in more detail and connect it to the evolution of the structure of the corresponding modes.

C. Eigenvector mixing and avoided crossings
Our conjecture is that inter-area oscillations directly originate from the set of zero-modes just described, corresponding to an appropriately chosen network partition. One key point is to show that, at least for few of the lowest-lying hybridized zero-modes diagonalizing the projection of L 1 onto the degenerate zero-subspace D, the linear combination is only weakly sensitive to the connection parameter ε ∈ [0, 1], i.e., well beyond the expected validity range of both perturbation theory and the standard theory of inter-area oscillations. To qualitatively understand why that is so, we first recall that, unless some underlying symmetry is present, the eigenvalues of a parameter-dependent matrix such as L(ε) generically do not cross as ε is varied. This resistance of eigenvalues to crossings is illustrated in Fig. 3 which shows that eigenvalues of L(ε) for a seven-area partition of the PanTaGruEl model exhibit avoided crossings [27] -they can get very close to one another but generically do not cross. To further demonstrate this, zoom-ins on three characteristic avoided crossings are shown in the top right panels of Fig. 3. Second, we recall von Neumann and Wigner's argument that an eigenvector of a parameter-dependent matrix remains essentially the same as this parameter is varied, as long as its eigenvalue stays away from avoided crossings [28]. We conclude then that, if an eigenvalue L(ε) does not go through an avoided crossing as ε increases from 0 to 1, the structure of the corresponding eigenvector does not change much This argument is corroborated by the data in Fig. 3. The change in eigenvector structure can be measured by the overlap η α (ε) = |u α (0)·u α (ε)| of an eigenvector at ε = 0 and at finite ε. Fig. 3 clearly shows that the lowest two hybridized modes (green and red curves) do not change their structure. The effect of an avoided crossing on eigenmode structure is clearly illustrated in Fig. 3. In the top panel, one sees that the seventh hybridized (blue curve) and first non-hybridized eigenmodes get close to each other at around ε 0.18 (indicated by the red circle), but do not cross. Simultaneously, a drop in η for the seventh hybridized eigenmode is observed in the bottom panel at the same value of ε. This reflects a scrambling of the eigenmode due to the avoided crossing. The same behavior is observed for the fourth (violet) and fifth (beige) hybridized modes, coincidentally at about the same value of ε.
The occurrence of avoided crossings as ε increases signals the onset of eigenvector mixing, beyond which the eigenvectors of L(0) are no longer representative of the eigenvectors of L(ε). Before we discuss such occurrences, we first derive a general criterion for the breakdown of perturbation theory for eigenvalues.

D. Perturbation Series Convergence Criteria
According to d'Alembert's ratio criterion the convergence radius r of a series f (x) = ∞ n=0 c n x n is given by r = lim n→∞ | cn cn+1 |. For the perturbation expansion of Eqs. (11) to converge up to ε → 1, the corrections need to be smaller with each order. For the eigenvalue expansion this translates into We approximate the expression on the left-hand side, assuming that u value for α ∈ D and for all β / ∈ D. Under this assumption, Eq. (12) reads We can get rid of the sum over β using from which we obtain, with a little bit of algebra, This expression is helpful, because it no longer contains a sum over β / ∈ D and expresses λ (2) α only as a function of λ (1) α and the expectation value of the squared interaction Laplacian over the eigenvector u (0) α . Note that the latter, despite fulfilling (15) are not eigenvectors of L 1 , so that u α L 2 The second order correction finally becomes where we introduced the generalized Kirchhoff indices of the disconnected subgraphs [29], From this analysis we find the following criterion for convergence The criterion combines an overall network characteristic with a mode-specific characteristic. To be satisfied, Eq. (27) requires that either (i) 1, or both. The first condition, that the weighted sum of subgraph Kirchhoff indices is small, requires that intraconnections are strong within subgraphs. This is so, because the Kirchhoff index is the average resistance distance in the graph. This first condition is consistent with recent works which find that the coherence of networks increases with the first non-zero eigenvalue of the network Laplacian, even for higher-order generator models [30]. The second condition requires that the inter-area connection strength felt by the α th hybridized mode is small. In this sense, Eq. (27) is qualitatively similar to the conditions for validity of singular perturbation theory and time-scale separation [5], [8], [9], [31], [32]. There is however a significant difference in that the condition of Eq. (27) applies to each hybridized zero-mode individually. In particular, perturbation theory may capture certain modes efficiently, while failing for others. This difference is of key importance, as we will see that for modes at the bottom edge of the spectrum -corresponding to the low frequency interarea oscillations -perturbation theory remains valid at higher ε, whereas the theory breaks down earlier for other modes higher up in the spectrum. The standard criteria for validity of singular perturbation theory and time-scale separation are global, therefore they implicitly request that the theory is applicable to all states. They are therefore too restrictive.
The criterion of Eq. (27) is helpful; however, it still needs to be computed numerically. It shows that the breakdown of perturbative approaches is mode-specific. Similar criteria can be found for higher orders of perturbation theory. Even though we cannot evaluate d'Alembert's criterion for n → ∞ we find that the first few orders already give a good approximation of the convergence.

E. Avoided Crossings
The criteria for validity of perturbation theory derived in the previous section are mode-dependent. They raise the important issue of determining which modes are best captured by perturbation theory for larger ε. To that end we recall what is sometimes referred to as the von Neumann-Wigner theorem [28], which states in our case that, as ε varies, pairs of eigenvalues of L(ε) = L 0 + εL 1 may undergo close encounters, however, they will generally avoid crossing each other's path. Von Neumann and Wigner further argued that it is at these avoided crossings that eigenvectors get mixed and that their structure changes fundamentally. Conversely, as long as an eigenmode is not undergoing any avoided crossing, then its structure does not change much. Therefore, predicting the first occurrence of an avoided crossing among the lowest eigenvectors is key to understand how far ε can grow, without altering the structure of an eigenmode. Stated otherwise, if the first few hybridized zero-modes do not undergo any avoided crossing until ε = 1, then the initial area partition at ε = 0 predicts the inter-area mode structure at ε = 1 with great precision. We analyze the situation for the case of two and of more than two areas.
1) Case of two areas: With two areas, there is one global zero-mode and one hybridized zero-mode. We want to derive a condition for the eigenvalue of the latter not to undergo an avoided crossing with the third eigenvalue, corresponding to the lowest non-zero-mode. To that end we first show that the first order correction λ (1) 2 to the hybridized zero-mode gives an upper bound to its true eigenvalue. Consider the first non-zero eigenvalue λ 2 (ε) at positionε andε + δε where 0 < δε 1. The slope at both points is given by the first order perturbation theory The eigenvector atε + δε can be expressed in terms of the eigenvectors atε using (13) Eq. (28b) then becomes The second term on the right-hand side of Eq. (30) is always negative, because the denominator of each term is negative while the numerator is positive. We therefore conclude where the lower bound is due to L 1 being Laplacian. Thus, λ 2 (ε) is concave and it is upper bounded by its first order perturbation theory correction at ε = 0. Note that Eq. (31) remains valid, regardless of the number of areas. Second, we derive a lower limit for the third smallest eigenvalue. Weyl's theorem [22] for the eigenvalues of the sum of two symmetric matrices states that With A = L 0 and B = εL 1 , this gives λ 3 (ε) ≥ λ 3 (ε = 0). Thus, a sufficient condition that there is no avoided crossing between λ 2 and λ 3 up to ε = 1 is This condition underestimates the validity of perturbation theory.
2) Case of more than two areas: When there are more than two areas, hybridized zero-modes may interact with one another via avoided crossings of their respective eigenvalues. The occurrence of these crossings is harder to predict than those between a single hybridized mode and the first non-zero mode treated in the previous paragraph. Here we focus on the occurrence of an avoided crossing between the first and second hybridized zero-mode. Eq. (31) remains valid and λ 2 is a concave and monotonously increasing function of ε. This behavior is captured only at small enough values of ε by a truncated perturbative series. This is not too big a restriction, since one expects avoided crossings to occur at low values of ε, because it is there that the largest changes to the topology of the network occur -going from unconnected to connected. Here we restrict our discussions to series truncated at (and including) second order in ε and construct conditions under which there is no avoided crossing between λ 2 and λ 3 before (2) 2 , where the truncated series reaches its maximum.
There are then two separate conditions under which there is no avoided crossing between λ 2 and λ 3 . The first one is when λ 3 , because then second order corrections increase the already increasing distance between λ 2 and λ 3 in first order perturbation theory. The second one is when the second order corrections are not sufficient to induce a crossing between the series for λ 2 and λ 3 truncated at second order before ε max . These two conditions read While the argument above leads to a = 1, we found numerically that a value of a = 1.5 gives better predictions.

IV. NUMERICAL VALIDATION
We validate the perturbation theory presented above by numerical investigations on three different networks: (i) a synthetic two area network, (ii) the IEEE RTS 96 test system, and (iii) the PanTaGruEl model of the synchronous grid of continental Europe.

A. Synthetic Two Area Network
We generate two Erdős-Rényi graphs, each with 50 nodes and a connection probability of p = 0.1, resulting in each node being connected to a bit less than 5 other nodes on average. All lines in these graphs have the same capacity, B ij ≡ 1, and we next connect them via (i) 5 and (ii) 25 lines, each with the same capacity B ij ≡ ε ∈ [0, 1]. The networks are shown in Fig. 4. The additional connections change the number of lines from 228 to 233 in the first case and to 253 in the second case. In both instances, two areas are still clearly defined, yet in the second case, we will see that the occurrence of an avoided crossing as ε increases totally mixes the structure of the single hybridized zero mode in this two-area set-up.
In the first case, we introduce five random connections between the two areas. We find that the left-hand side in Eq. (27) is significantly smaller than one, so that perturbation theory should be valid and the slowest inter-area mode should reflect the structure of the network. This analysis is confirmed by calculating the actual perturbational corrections up to third order (not discussed above, for details, see Ref. [22]) and noticing that with each order the approximation converges to larger ε. We furthermore check that Eq. (33) is satisfied, therefore we expect that the first eigenmode is well captured. This is the case with η(ε = 1) = 0.99 indicating that the eigenvector barely changes with increasing ε. The top right panel in Fig. 4 shows that in this case, λ 2 (ε) undergoes no avoided crossing.
In the second case, twenty-five random connections are introduced between the two areas, so that every other node has a connection to the other area on average. The bottom right panel in Fig. 4 shows that, again in this case, perturbation theory approximates the eigenvalue well even for ε → 1. However, this time we find that Eq. (33) is not met. We therefore expect the lowest eigenmode to undergo an avoided crossing, and this is confirmed in the bottom right panel of Fig. 4, where an avoided crossing is visible at ε ≈ 0.7.

B. IEEE RTS 96 test system
The IEEE RTS 96 test system consists of 73 nodes divided into three well-defined areas as shown in Fig. 5 [33]. We consider two different initial aggregations, first along the obvious area boundaries, second along three lines cutting across each area.
In the first case, we find that the slow eigenvalues are well captured by perturbation theory, and that moreover there is no avoided crossing affecting the two lowest non-zero eigenvalues corresponding to the hybridized zero-modes. Accordingly, the corresponding modes retain the structure of the initial aggregation.
In the second case, the chosen initial aggregation results in large perturbative corrections and eventually to the lowest eigenvalues undergoing an avoided crossing, as predicted by Eq. (34). The avoided crossing is shown in the bottom right panel of Fig. 5. This shows that an incorrect initial aggregation leads to avoided crossings, which are the mechanism for mode mixing. Numerical investigations show that indeed the structure of both eigenvectors arising from the degenerate subspace changes almost completely (η(ε = 1) < 0.06).

C. PanTaGruEl
The PanTaGruEl model of the synchronous grid of continental Europe consists of 3809 nodes connected by 4944 power lines [11], [21]. With the dispatch used for this paper, there are 468 generators. To illustrate the validity of the theory presented above, we used the standard aggregation algorithm of Ref. [31]. We found that, to capture the interarea oscillations, an aggregation into seven areas works well, and that the number and structure of these modes does not change as the number of areas increases.
PanTaGruEl is a strongly connected network with no obvious area separation (except perhaps the Iberian Peninsula) and, not surprisingly, the convergence criterion (27) is by far not met. Yet, the slowest eigenvectors are not much affected by the inter-area connections. This is predicted by the criteria , and a counterintuitive aggregation (red, blue and green nodes). b) ε-dependence of the three lowest eigenvalues for the correct aggregation. There is no avoided crossing and the eigenmode preserve their structure all the way to ε = 1. c) Color-coded structure of the lowest mode of the Laplacian for ε = 1. Its structure is already well predicted by our first-order perturbation theory with the correct aggregation, giving an overlap η 2 (ε) = |u 2 (ε = 0) · u 2 (ε = 1)| 0.96. d) ε-dependence of the three lowest eigenvalues for the counterintuitive aggregation. There is an avoided crossing between the second and the third eigenvalues and their structure is significantly changed well before ε = 1. The second order corrections are shown as dashed lines to visualize Eq. (34).
of Eqs. (34), which are met, and corroborated by numerical data. First, Fig. 6 shows the lowest nonzero eigenvector of the network Laplacian at ε = 0.1 and ε = 1. Clearly the general structure remains the same, regardless of the interarea coupling strength. This is corroborated by the overlap data shown in Fig. 3. Interesting is the behavior of the fourth and fifth eigen-  vectors which undergo an avoided crossing shown in Fig. 3 at ε 0.2. Their behavior illustrates the eigenvector mixing process discussed above, as after the avoided crossing, ε > 0.2, the actual eigenvectors are given by u 4,5 ≈ (ũ 4 ±ũ 5 )/ √ 2, whereũ i denotes the eigenvector before the avoided crossing. This results in overlaps η 4,5 (ε = 1) ≈ 0.45 at ε = 1 for both eigenvectors. The seventh eigenvector also undergoes an avoided crossing at ε 0.2, which mixes it with the high-lying eigenvectors. From Fig. 3 we see that this is accompanied by an abrupt decrease of η 7 . These phenomena nicely illustrate the direct connection between avoided crossings and changes in the eigenvectors structure. Figs. 1 and 6 show that the lowest hybridized zero-mode essentially resides in the Iberian Peninsula and in the Balkans. These are two of the seven areas in our aggregation. We found that oscillations between these two areas are triggered by a perturbation in either one. This is shown in Fig. 1 in the case of a noisy power injection. The chosen noise is an Ornstein-Uhlenbeck noise with a correlation time much larger than the oscillation frequency, and we found that the latter is close to the oscillation frequencies of the two lowest hybridized zeromodes and not related to any perturbation time scale. These two modes are excited by the perturbation, and the addition of their components evidently resolves the two areas as shown in Fig. 1.
As another example, we finally investigate the reaction of the system to a 900 MW power loss. Fig. 7 shows the response of the Balkan area and the Iberian Peninsula to a fault in the opposing area. As before, all nodes within one area respond coherently -with the same frequency and phase -to a fault in the opposing area, with only moderate variations of the response amplitude. The Fourier transform of the oscillation response in the Iberian Peninsula is further shown in the bottom panel. Superimposed on it are the locations of the eigenfrequencies of the A-matrix, and the Fourier spectrum indicates that only the first two eigenmodes are excited, with a broadening originating from the damping in Eq. (4). This confirms our above claim that inter-area oscillations in the PanTaGruEl network are mainly carried by the two lowest modes of the system.

V. CONCLUSION
The standard approach to slow coherency successfully predicts the structure of slowly oscillating, long-wavelength interarea modes, even in strongly connected power networks that lie outside its range of validity. The theory presented above puts this theory on solid grounds even in such well connected networks. Our line of reasoning goes as follows.
First, we recalled that in homogeneous systems with constant inertia and damping, small-signal oscillations in the swing equations (1) are carried by eigenmodes of the A-matrix of Eq. (6). The structure of these modes is solely determined by the eigenmodes of the Laplacian matrix of the network [20].
Second, when dealing with slow coherency/inter-area oscillations, the focus is on the slowly oscillating modes. A recent extension of Courant's nodal domain theorem [13] shows that the slow modes have large nodal domains [14]. Because of that, these slow modes are only poorly resolving inhomogeneities in inertia and damping and are therefore much less sensitive to them. Perturbation theory makes this quantitative and shows that, while in the presence of significant inhomogeneities, most eigenmodes of the A-matrix acquire a structure not necessarily captured by those of the network Laplacian L. The long-wavelength modes of A retain the structure of the slow modes of L.
Third, using perturbation theory, we showed that an appropriate choice of disaggregation of the network into initially disconnected areas captures the slow modes from the hybridization of the zero modes of each area Laplacian, even when the inter-area couplings are restored. This clarifies the origin of the slow coherency, inter-area modes, and in particular explains how they often are quasi-homogeneous over large areas -they originate from area-Laplacian zero-modes that are exactly constant there.
The aim of our theory was not to identify the optimal initial aggregation. We found numerically that standard numerical aggregation algorithms do a good job at predicting that. Our theory fills an important theoretical gap, in that it explains (i) the origin of the slow inter-area modes and (ii) why the standard approach to inter-area oscillations still works well outside its range of validity. Doing so, we closed a number of loopholes in the theory of inter-area oscillations in transmission power networks.