Propagation of Disturbances in AC Electricity Grids

The energy transition towards high shares of renewable energy will affect the stability of electricity grids in many ways. Here, we aim to study its impact on propagation of disturbances by solving nonlinear swing equations describing coupled rotating masses of synchronous generators and motors on different grid topologies. We consider a tree, a square grid and as a real grid topology, the german transmission grid. We identify ranges of parameters with different transient dynamics: the disturbance decays exponentially in time, superimposed by oscillations with the fast decay rate of a single node, or with a smaller decay rate without oscillations. Most remarkably, as the grid inertia is lowered, nodes may become correlated, slowing down the propagation from ballistic to diffusive motion, decaying with a power law in time. Applying linear response theory we show that tree grids have a spectral gap leading to exponential relaxation as protected by topology and independent on grid size. Meshed grids are found to have a spectral gap which decreases with increasing grid size, leading to slow power law relaxation and collective diffusive propagation of disturbances. We conclude by discussing consequences if no measures are undertaken to preserve the grid inertia in the energy transition.

In order to cover the increasing human energy demand by renewable energy resources and to ensure that this energy will be available wherever and whenever it is needed, more efficient energy transport and storage technologies need to be developed. The fluctuations in generated power by wind turbines and solar cells -both in time and geographically -demand to explore new strategies to store energy on all time scales and to distribute the power in the grid smartly. At the same time, the spreading of critical disturbances throughout the grid has to be prevented to ensure the stability of the entire grid. Renewable energy resources fluctuate strongly in time on time scales as small as seconds. Moreover the inverter-connected wind turbines and solar cells provide no inertia 1 . This is in contrast to conventional generators, whose rotating masses hold inertia and thereby momentary power reserve available for the grid, which makes the grid resilient and prevents strong fluctuations of the grid frequency on time scales of several seconds 2,3 . As the inertia in the grid keeps decreasing with higher share of renewables, the grid is responding on shorter time scales to disturbances. It is therefore essential to understand the impact of this development on the stability of electricity grids. In this article, we aim to find out if and how the relaxation and propagation of disturbances in AC grids is modified when the grid inertia from the rotating masses of generators is decreasing. We focus on the spreading of weak disturbances, which do not cause by themselves a failure of transmission lines or generators.
The dynamic interaction and response of generators and consumers is studied modeling the grid as a network of nonlinear oscillators [2][3][4][5][6][7] . These nonlinear dynamic power balance equations, so called swing equations, describe the dynamics of coupled rotating masses by a system of coupled differential equations for local rotor angles φ i , where i denote the grid nodes. We aim to analyze the propagation of disturbances in AC grids on short time scales up to several seconds. Therefore, we do not consider control measures, which typically set in on longer time scales 3 . We expect that a better understanding of this short time transient dynamics will then contribute to the optimization of primary and secondary control measures [1][2][3] . Furthermore, we expect that the understanding of the spreading of small disturbances will be important for the design of improved power system stabilizers (PSS). Thereby it can contribute to prevent the occurrence of larger disturbances, whose spreading could result in overload and failure of transmission lines, and cause cascading failures. We note that the dynamics of cascading failures involves additional mechanisms, such as the overheating of transmission lines. While most studies of cascading failures are based on stationary power balance calculations or simplified topological flow models such as the messenger model [8][9][10][11][12][13][14] and our own work based on nonlinear power flow calculations 15 , only few model their dynamics 16,17 , introducing for example a time scale for overheating of transmission lines.
In this article, we solve the nonlinear swing equations numerically and explore how a local perturbation propagates throughout the grid. The origin of disturbances can be fluctuations in generating power or sudden changes of transmission line capacitance. We analyze these results, employing analytical results obtained from a linear response theory, mapping the swing equations on discrete linear wave equations for small perturbations 18 . Depending on the geographical distribution of power, power transmission capacity and grid topology we find that the disturbance may either decay exponentially in time with the decay rate of a single oscillator Γ 0 , or exponentially with a smaller decay rate Γ < Γ 0 , or, even more slowly, decaying with a power law in time. Such a slow power law decay is found to arise together with slow, diffusive propagation 18 .

Model Description
Grid Topologies. Aiming at a systematic approach, we consider three different grid topologies. The Cayley tree, Fig. 1(a), resembles distribution grids which are operated in a tree-like fashion to pinpoint and repair failures more easily. Branches grow outward from a central node with branching number b, forming l branching levels. Tree grids are characterized by the distance between neighbored nodes a, the total number of nodes N, branching number b and level l. The degree d i , the number of links connecting node i to any other node is an important network characteristic. For Cayley trees it is d = b + 1 except at outer edge nodes where d = 1. The square grid, Fig. 1(b), is a meshed grid, used as basic model for transmission grids with their strict redundancy demand to guarantee continuing operation when a single line fails (n-1 criterion), characterized by distance a, linear grid size L and number of nodes N = (L/a) 2 . Their degree is d = 4, except at edges (d = 3) and corners (d = 2). Thirdly, we use the open-source SciGRID dataset 19 of the German transmission grid, Fig. 1(c), as a real-world example of a highly meshed grid. Excluding island nodes, the largest connected network of the three highest voltage levels, 400 kV, 380 kV, 220 kV and some 110 kV lines has N = 502 nodes and 673 links. Its degree d i has a wide distribution, Fig. 2, with average degree 〈d i 〉 = 2.7 and typical degree = = . .
Excluding singly connected nodes (which are mostly artifacts since the data set does not include the transnational European grid) we get an average degree 〈d i 〉 = 3.5.
Dynamic AC Grid Model. AC transmission grids are three-phase systems which are typically loaded and operated symmetrically, so that the power flow balance equations depend on a single phase, only. Neglecting Ohmic losses along the lines, which are small in high voltage transmission grids, we assume purely inductive transmission lines with combined inductance L ij between nodes i and j with power capacity K ij = |V i ||V j |/(ωL ij ) where V i , V j are the voltages at nodes i and j. ω = 2πf with grid frequency f = 50 Hz. We assume constant voltage amplitude V throughout the grid, V i = V exp(ıΦ i ). Since voltage amplitudes change typically on larger time scales than phases, we will focus on the dynamics of phases Φ i . Fixing the voltage V there are no dynamic terms in the reactive power balance equation (they appear in higher order when voltage dynamics in addition to the phase dynamics is considered 3 ), so that we need to consider active power balance equations, only. Since our main goal is to study the influence of grid topology and inertia on the phase dynamics we assume equal inductance L and power capacity K = V 2 /(ωL), for all links, yielding the power capacity between nodes i and j as K ij = KA ij , where A ij is the grid adjacency matrix. Thereby, the stationary active power flow balance equations are obtained from Kirchhoff 's laws as 20 However, as loads and generated power vary in time, that power balance may become violated and the nodal phases become dynamic. Denoting the solution of the stationary balance equation, equation (1) by the phase shift θ , i 0 we can write where α i (t) are the dynamic phase shifts. The grid nodes are either connected to synchronous generators with inertia J i or to loads which can be motors with a finite inertia, being modeled as synchronous motors 3,21 . The electric power P i , is either positive (generator) or negative (motor). In order to be able to study the dependence on system parameters and topology we take homogenous parameters P i = s i P, with s i ∈ {+, −} and J i = J. The phase dynamics is governed by the balance of changes in kinetic energy, energy dissipation and electric power exchange with adjacent grid nodes, yielding the swing equations [4][5][6]21 .
where γ is a damping constant, F ij = K ij sin(Φ i − Φ j ) the power flow in the transmission line between nodes i and j.
If the nodes would not be coupled by the transmission lines, the phase at each node would decay exponentially fast with the local relaxation time τ = J/γ, which increases with inertia J and decreases with damping parameter γ. Therefore, when studying the temporal and spatial dependence of α i , it is convenient to scale the time with relaxation time τ. Inserting equation (2)

Results
Transient Dynamics. In order to study the transient behavior of AC grids perturbed by local disturbances, we solve the nonlinear swing equation (4) on the different grid topologies of Fig. 1 as function of the set of parameters (τ,Π K , Π P ).

Stationary solution.
Before any perturbation is applied, we calculate the stationary state θ i 0 at every node i in the grid. This is accomplished by first obtaining a solution of equation (1) for small phase differences linearizing where P and θ 0 are vectors, whose i -th component is the power and stationary phase at node i, respectively. H has at least one Eigenvalue zero. Therefore, we need the pseudo inverse H + , yielding θ 0 = H + ⋅ P. We use this solution as initial condition for a numerical root solver to find the solution of the nonlinear equation, equation (1). This way, the numerical accuracy of the stationary solution is maximized to make sure that we use as initial condition for the swing equations the stationary state.  (1), we insert them into swing equation (4). Next, we solve these equations when the AC grid is perturbed by a local disturbance as outlined in detail in Supplementary I. As disturbance of the stationary state we increase the power for a short time interval 0 ≤ t ≤ Δt pert at the grid node marked with 'x' Fig. 1. We choose as small perturbing power one per mille of the initial generator power P. The resulting transient behavior of phase deviations α(t) is shown in Fig. 3 for the Cayley tree grid of Fig. 1(a) as function of rescaled time t * = t/τ for parameters (Π K = 10, σ = Π P /Π K = 0.08). We define the distance between any two nodes ⁎ r ij as the number of lines in the shortest path connecting them 22 . We see that the phase at the disturbed node, r * = 0 is perturbed first, reaching a maximum after a delay time, decaying then in an oscillatory manner. Phases of nodes further away from the origin of the disturbance are perturbed later and reach smaller amplitudes. We analyze this temporal and spatial propagation of the perturbation quantitatively below. Since we are interested in the propagation of small disturbances which do not destabilize the system, we review next the conditions for stability in order to make sure that we choose the size of the disturbance accordingly.
Stability. Without disturbance, the criterion for stable (allowed) and unstable (forbidden) parameters Π P , Π K is the existence of a non-complex solution to stationary state equation (1). Thus, the ratio σ = Π P /Π K determines whether parameters are allowed. A critical value σ c exists which may depend on grid topology and power distribution. If there are no clusters of consumers or clusters of generators, the critical value at node i is given by σ ci = d i , where d i is the node degree. Thus, the critical value, below which the whole grid is stable, is given by σ c = min(d i ). For a general distribution of P i there can be clusters of generators or consumers and the critical value σ c is determined by the size and form of these clusters. If a cluster of generators has total power P C = ∑ cluster P i , with effective degree d C , as obtained by counting the number of consumers which are directly coupled to that cluster, the critical value is given by σ c = min(P d C /P C ). As outlined in the article in Supplementary section II the damping term and therefore the parameter τ = J/γ determines the size of the basin of attraction, if a stable fixed point exists, that is for σ < σ c . Outside of the basin of attraction the system is attracted to unstable fixed points, so called limit cycles. For large damping there can be a regime where there is no coexistence with limit cycles and the whole phase space is stable, see f.e. refs 23,24 for a review. Therefore, depending on the magnitude α of the perturbation it can destabilize the grid already at σ * (α) < σ c . In the Supplementary section II. we derive a typical upper limit for the size of the perturbation α before it kicks the system out of the stable region and find that the disturbance indeed destabilizes the grid already at values σ * (α) < σ c . Only in the limit, when the perturbation amplitude is vanishing we recover σ * (α → 0) = σ c .
Classification of Transient Dynamics: Parametric Phase Diagrams. In order to analyze the transients we calculate absolute values of the power flow change in the transmission line between nodes i and j, as averaged over all lines at a distance r from the disturbance and divided by its maximum value, In Fig. 4 we show examples for transient dynamics Δf (r * , t * ) in three grid topologies for different values Π K , σ. We next varied parameters Π K and σ in small steps. Identifying all parameters with unstable solutions, we find unstable parameter regions σ σ > , c shown in red for the tree grid in Fig. 5(a), in (b) for the square grid with random arrangement of consumers and generators (the result for a periodic arrangement is shown in Supplementary  Fig. 3  random arrangement and σ c = 0.33 for the German transmission grid. In stable regions we identify three qualitatively different transient behaviors: The perturbation may decay exponentially fast (FE) with local relaxation rate Γ 0 superimposed by oscillations, as seen in Fig. 4(a) in the tree grid. All parameter sets showing FE behavior are plotted in Fig. 5(a-c) as green circles. More examples for fast exponential transients are shown in Supplementary  Fig. 2(d-f). Secondly, we observe exponential decay with a smaller relaxation rate Γ < Γ 0 for a large interval of time τ −  t t 0 as seen in Fig. 4(b) for the tree grid, d) for the square grid and f) for the German grid for the parameter sets shown in the phase diagrams Fig. 5(a-c) as yellow squares. Green shading denotes areas where fast relaxation (FE) is expected to occur, as limited by the dashed lines in Fig. 5(a,b), corresponding to the analytical result, equation (16) for the tree grid and to equation (11) for the square grids, respectively. For the German grid, we indicate the boundary of that region in Fig. 5(c) by a dotted line, according to numerical results. The yellow shaded areas in Fig. 5(a,b) show where slow relaxation (SE) is expected to occur, in good agreement with the analytical line. The slight inconsistency is in a region where we observe only small deviations of Γ from Γ 0 well within the estimated error bars of the numerical results. We plot the parametric dependence of the relaxation rate Γ(σ) in Fig. 6 for b = 3 Cayley tree in units of Γ 0 , as obtained by fitting the transients with an exponential decay (error bars denote the range of Γ which provided good fitting to an exponential decay). Thirdly, we observe an even slower decay with a power law in time in square grids and in the German grid, as seen in Fig. 4(c) for a square grid and Fig. 4(e) for the German grid. That power is found to be close to 2, which is in very good agreement with the diffusive behavior of the change in power flow Δf(t) obtained for a square grid by an analytical derivation 18 , as reviewed in the next section, equation (13). The numerical results indicate that such diffusive behavior may occur also in other meshed grids, as the example of the German grid at small values of the parameter Π K < 2 shows, Fig. 4(e). Thus, as the inertia J and thereby Π K is lowered, more nodes become correlated and the spreading of a disturbance is slowed down to a collective diffusive spreading for times τ > t .
Spatial Propagation of the Disturbance. Next, we examine the spatial propagation of disturbances. If it propagates ballistically, it moves with velocity v, so that it reaches a node at distance r after time t according to r = v(t − t 0 ), where t 0 is the time when it occurs first. Diffusion corresponds to areal spreading, so that a node at distance r is reached according to r 2 = 4D(t − t 0 ) with diffusion constant D. Localization occurs when the disturbance does not spread beyond a certain localization length ξ, never reaching nodes at distances ξ > . r Let us start by analyzing the propagation in square grids with L = 22 and random arrangement of generators and motors. In Fig. 7 we show numerical results for σ = 0.1 for the arrival time t * = t/τ (after initial disturbance δΠ P = 0.001Π K occurs at time t 0 = 0s), which the power disturbance needs to reach nodes at geometrical distance r * = r/a, where = + r r r ( ) x y 2 21/2 (results for distances r, defined as the minimal number of links connecting two nodes, are shown in the Supplementary). We define the arrival time as the time, when the phase deviation α i exceeds a threshold α th = 10 −6 δΠ P at a node at distance r. In Fig. 7(a) we plot results for Π K = 10 5 together with t * = cr * , where c = 0.0034 is fitted (red line). According to the analytical result it propagates for Π > Π K K s ballistically with the velocity given by equation (12), as derived in the next section. We plot t * = r * /(vτ/a), as the dotted red line. Thus, we can confirm quantitative agreement with ballistic motion. Here, Π K s is given for the square grid by equation (11), as plotted in the phase diagram Fig. 5, the dashed line. b) for smaller inertia, Π K = 10, we observe a crossover behavior where the data in Fig. 7(b) fits neither ballistic nor diffusive propagation well. For Π K = 0.1, we find in Fig. 7(c), that the arrival time agrees with diffusive motion, t * = cr *2 where we fitted c = 1.1083 (pink line), while the ballistic formula (straight lines) do not fit the data. Using the analytical result for the arrival time t * = r *2 /(4D)f th (r), where f th (r) is a logarithmic correction, which depends on threshold value α th , as derived in Supplementary III, we find that it provides a good lower bound to the data (pink dashed line).
In Fig. 8 we show for the german transmission grid the time t * = t/τ which a power disturbance, initially at the node marked by a cross in Fig. 1, needs to reach a node at distance r * = r/a. For σ = 0.1 and Π K = 10 5 we fit the numerical results with the ballistic formula t * = cr * with slope c = 0.0062 (red line), and compare it with the analytical value for velocity v, equation (12) (dashed red line). We see that the analytical formula provides a lower bound to the data, yielding evidence for ballistic motion for some nodes, until the node distance r = 9, which corresponds to the shortest path length from the source of the disturbance to the boundary of the german grid. The disturbance is seen to reach other nodes later, indicating anisotropy in the german grid. For (b) Π K = 10 we show the fitted curves for ballistic (red line) and diffusive motion (pink line) and analytical curves (red and pink dashed lines). While diffusive spreading is fitting some nodes, the disturbance needs more time to reach other nodes, which is another consequence of anisotropy. For (c) Π K = 0.1, corresponding to low inertia in the grid, we show the analytical curves for ballistic (red dashed line) and diffusive motion (pink dashed line). The scatter of the results for different nodes at distance r is too large for a meaningful fit. The analytical formula for diffusive spreading is seen to provide a good lower bound for the data. Thus, we find strong indications in the german  grid that for low inertia the collective dynamics of nodes results in diffusive spreading of the disturbance. The spreading is more strongly delayed for some nodes, and there are indications of localization of disturbances in certain directions since some nodes do not become excited above the threshold within the observation time. We also note that for some nodes the relaxation to a stationary state did not take place before the signal started, see the Supplementary material for more details.
On the Cayley tree grid the time t * = t/τ the power disturbance needs to reach a node at distance r * = r/a for power ratio σ = 0.1 is shown in Fig. 9. For a) Π K = 10 5 we fitted the numerical results with the ballistic formula t * = r * /(vτ/a) (pink line), and the quadratic formula (dashed pink line). As shown in the next section, the quadratic dependence t * = cr *2 , where c is a constant, is due to the quadratic dispersion of the Eigen modes of the generalized Laplacian on the tree grid. It fits for b) Π K = 10 the numerical data even better. Note that this diffusion like dependence occurs at values of Π K , where we observed fast exponential decay. This is explained in the next section. Note also that for both Π K = 10 5 and Π K = 10 the spreading is isotropic in all directions, as the comparison between node resolved (full symbols) and same distance r node averaged (white symbols) shows. As the inertia is lowered, there appears a strong anisotropy at large distance, however, as the results for Π K = 0.1 show. Thus, we find for small inertia anisotropic spreading, while the disturbance still decays exponentially in time. This is explained in the next section by a topologically protected spectral gap in treelike grids.

Linear Response theory and Spectral Analysis of Disturbances: Analytical results. For small
perturbations, corresponding to the condition σ < σ * (α), we analyze the swing equation (4) by expanding it in the perturbation α i − α j . This yields the linear wave equations on the grid 18 , ij Kij i j 0 0 depending both on grid topology and on initial distribution of power P i through stationary phases θ . i 0 These coupling amplitudes define the generalized Laplace operator Λ with Λ ij = −t ij and Λ ii = ∑ i t ii , which is related to the stability matrix in linear stability analysis 21,25 as used in small signal stability analysis 26,27 . We expand the phase deviation α i (t) and the disturbance δΠ i (t) = JδP i (t)/(γ 2 ω) in a generalized Fourier series in terms of the Eigenvectors φ n of Laplacian matrix Λ, as defined by Λφ n = Λ n φ n , where Λ n is its Eigenvalue 18,21,25,28 , see Supplementary III. For a local perturbation at a site j, lasting a short time interval and Eigenvectors φ n of the Laplacian, which can be obtained numerically by exact diagonalization for arbitrary grids. For particular grids we obtained analytical solutions as summarized in the following.
Periodic Square Lattice. For square lattices where P i = ±P are arranged periodically, an analytical solution of the stationary power balance equation is obtained with plain wave Ansatz φ = c e qi q iqr i with wave vector q. The Eigenfrequencies ε q n are 18 , For finite size L, the wave vectors are quantized, q x,yn = n x,y π/L, with n x,y = 0, ±1, ... . The resulting dispersion of the Eigenfrequency ε q is plotted in Fig. 10 as function of discrete wave numbers q n (blue dots). We observe a spectral gap to the first excitation energy ε = Δ q L 1 as indicated by the dotted line in Fig. 10,   Figure 9. The time t * = t/τ the power disturbance needs to reach a node at distance r * = r/a for exemplary sets of parameters in a Cayley tree grid. The initial disturbance δΠ = 0.001Π K occurs at t 0 = 0s.
SCIEnTIfIC REPORTS | (2018) 8:6459 | DOI:10.1038/s41598-018-24685-5 L K c 1/2 2 2 1/4 0 decreasing with size L. Applying linear response theory for small disturbances, inserting τ ε Λ = n q 2 2 n into the Fourier expansion, equation (8), we get the phase deviation α(t) for all times t. The condition that slow modes with small relaxation rate Γ < Γ 0 . Appear is, that the spectral gap Δ L is smaller than the local relaxation rate, Γ 0 . This yields the parametric condition Π < Π L ( ), depends on grid size L and power ratio σ as is fulfilled, Eigen modes with small wave number q have purely imaginary complex Eigenfrequency Ω q which results in slow decay of the phase deviations without any oscillations. Summing over all slow modes in the spectral representation of α i (t), equation (8) we find that a perturbation at node j propagates for times τ > t and distances exceeding the mean free path l = vτ diffusively with diffusion constant D = v 2 τ, see Supplementary III for the derivation. Diffusion causes slow power law relaxation of the disturbance at the initial site and an initial increase, followed by a power law decay at other sites. The resulting power law relaxation of the change in transmitted power between nodes k and l is 18 kl kl Thus, we find that the change in transmitted electric power at the site of perturbation r j decays with a power law in time with power 2 in excellent agreement with the numerical results for the periodic square grid, Fig. 4(c). The position of the maximum of the disturbance spreads according to = r D t 4 max 2 with time. Note that r is here defined as the geometrical distance = + r r r ( ) given by equation (11), where σ c is the critical value for that distribution of power P i . Diffusion occurs with an accordingly modified diffusion constant D. We plot Π L ( ), K s dashed line in Fig. 5(c), together with numerical results where we use the numerically obtained value for σ c .
Cayley Tree Grid. On a Cayley tree grid every inner node is connected to d = b + 1 other nodes, as shown in Fig. 1(a) for b = 3. For the Cayley tree with branching > b 1 symmetric eigenvectors were found in ref. 30 . With these Eigenvectors we obtain the Eigenfrequencies of the discrete wave equation 7 as can not be identified with a wave number since the phase of the Eigenvectors depends nonlinearly on qa. We plot ε q in Fig. 10 for a finite sized tree as function of discrete number q n . It is remarkable that Eigenfrequencies equation (14) have for > b 1 a finite gap Δ, independent on the number of nodes N, K c 1/2 2 2 1/4 0 indicated by the dashed line in Fig. 10. The condition that slow modes with Γ < Γ 0 , appear is Δ < Γ 0 , which yields the parametric condition Π < Π , If that condition is fulfilled, the perturbation decays for large times exponentially with a reduced relaxation rate Γ min = (1 − (1 − τ 2 Δ(z, Π K ) 2 ) 1/2 )Γ 0 which can be much smaller than the local relaxation rate Γ 0 . The spatial spreading is diffusive, s the perturbation decays fast exponentially with local rate Γ 0 = 1/τ. Due to the quadratic dispersion of ε q for a Cayley tree, at small values of q, see Fig. (10), the spatial spreading scales with Fig. 5(a) as the dashed line. For a tree grid with inhomogenous distribution of power P i , typically, the slowly decaying modes are expected to appear when Π < Π , K K s with Π K s given by equation (16), where σ c is the critical value for that distribution of power P i .

Discussion
In conclusion, we studied how the relaxation and propagation of disturbances in AC grids is modified when system parameters like the inertia in the grid are changed. To this end we solved the nonlinear dynamic power balance equations on three different grid topologies numerically and analyzed the results comparing them quantitatively with analytical insights obtained by linear response theory and spectral analysis. By a generalized Fourier expansion for the square grid and the Cayley tree grid we show that the long time transient behavior is governed by the spectral gap between the stationary state and the lowest Eigenmode of its grid. The Cayley tree grid has a finite spectral gap, which is protected by its topology and independent on grid size. Meshed grids, however are found to have a small spectral gap which is reduced strongly with increasing grid size, leading to slowed relaxation and collective diffusive propagation of disturbances. Analyzing the numerical results we confirm that, depending on inertia, geographical distribution of power, grid power capacity and topology, the disturbance may either decay exponentially in time with the decay rate of a single node, or exponentially with a smaller decay rate or, even more slowly, decaying with a power law in time, resulting in collective diffusive propagation. Let us discuss the relevance of our results for real grids like the german high voltage transmission grid, where the lines have a typical capacity of K ij = 10 GW 19 . Assuming that half of the nodes act as generators and the other half as consumers to meet Germany's peak power production of 83 GW 31 , we have |P i | = 300 MW. Typical conventional power plants of that rated power have inertia J = 10 4 kg m 2 and damping γω 2 = 0.10P. This yields Jω 3 = 310 GW, Π P = 1.03 ⋅ 10 5 and Π K = 3.44 ⋅ 10 6 . Taking the condition Π < Π L ( ),

K K S
with equation (16) as an estimate that meshed grids show diffusive behaviour, that condition becomes for currently existing transmission grids, > . L a 1856 Thus, diffusive propagation is unlikely to occur in present transmission grids even on the European scale. However, as conventional power plants become substituted by renewable energy resources the inertia in the grid is reduced substantially 1 . Thus, the dynamics of transmission grids will change. For very small inertia of J = 0.1 kg m 2 and otherwise same parameters, we find Π P = 1.03 and Π K = 34.45 so that the condition to observe diffusion is > . L a 5 87 and becomes relevant for transmission grids on the national scale. If no measures are undertaken to substitute the inertia of conventional power plants 32 , we conclude that the energy transition will change the transient dynamics drastically, disturbances will relax more slowly and spread by collective diffusive propagation, reducing the grid stability in meshed grids. Our results on tree grids suggest that a reduction of meshedness, may help to localize disturbances in the grid and yield faster relaxation, improving thereby the grid stability.
The linear response theory presented and applied here can be extended to larger perturbations by including nonlinear terms which introduce couplings between Eigenmodes of the Laplacian. This can be taken into account in various approximations, which will be considered in future research. Our results may also have important consequences for other systems modeled by a network of nonlinear oscillators, ranging from metabolic systems, to neuroscience and to supply chain networks. After submission of this manuscript, we learnt of a related work, which employs a linear response theory to calculate stability measures in 1st order Kuramoto models, corresponding to the swing model without inertia 33 .