Boosted fluctuation responses in power grids with active voltage dynamics

Secure electric energy supply and thus stable operation of power grids fundamentally relies on their capability to cope with fluctuations. Here, we study how active voltage dynamics impacts the collective response dynamics of networked power grids. We find that the systems driven by ongoing fluctuating inputs exhibit a bulk, a resonance, and a localized grid frequency response regime, as for static voltages. However, active voltage dynamics generically weakens the degree of localization in the grid, thereby intensifying and spatially extending the high-frequency responses. An analytic approximation scheme that takes into account shortest signal propagation paths among the voltage, phase angle and frequency variables result in an asymptotic lowest-order expansion that helps understanding the boosted high-frequency responses. These results moreover offer a generic tool to systematically investigate fluctuation response patterns in power grid models with and without active voltage dynamics.


Introduction
Climate change has been identified as 'a major concern for humankind' in the 2015 Paris agreement [1]. Its mitigation has thus been declared a major objective by the United Nations to reduce greenhouse gas emissions and thereby limit average temperature rise to 1.5 K. Achieving the objective of a significant reduction of greenhouse gas emissions, requires massive changes across all human activities, from agriculture and forestry, transportation, building infrastructures, heating and most large-scale industrial processes all the way to a reliable supply of electricity [2]. Such measures imply to decarbonize electric power generation and thus to replace large centralized power plants by power generation units based on renewable energies, such as solar-and wind power generators.
Renewable sources constitute smaller, more distributed and less predictable power inputs than traditional large plants, causing fluctuations in supplied power [3][4][5]. Moreover, increased digitization and changes of behavioral patterns at consumers of electric power cause more fluctuating power demands [6,7].
Fluctuations in the collective dynamics of power grids specifically originate from fluctuating power inputs due to rapidly changing local weather conditions and from decreasing the fraction of mechanical inertia in the generators [8,9] that are capable of dampening fluctuations. Fluctuating power demand in parallel to increasing shares of renewable energy sources thus induces stronger and stronger fluctuations acting upon operational steady states in electric energy transmission systems. To plan and operate future-compliant grids, we thus need to thoroughly understand their distributed response dynamics to external fluctuations.
In this article, we study the impact of active dynamics of voltage amplitudes in power grid models onto their collective response dynamics to fluctuating signals. We extend research on the standard second-order network models of power grids that incorporate rotor angle and frequency dynamics [10] to third-order models, the lowest order models that also capture voltage amplitude dynamics. Direct numerical simulations of response dynamics of these two model classes indicate that the collective responses of both model classes are roughly similar and exhibit three different characteristic regimes for rotor angle and frequency responses-homogeneous bulk responses at low frequencies, heterogeneous resonant responses at intermediate, and localized responses at high frequencies.
Intriguingly, however, active voltage dynamics boosts the high-frequency responses such that responses stay less localized and may have substantially larger amplitudes at specific response nodes in the network. Analytically deriving a network linear response theory for the third-order model, we reveal the origin of this boost. The dynamic voltage variables offer an effective shortcut along which frequency and phase angle variables of one node dynamically influence other nodes in the grid.
Based on the results of linear response theory, we analytically demonstrate that and how strongly these shortcuts increase the effective influence range of high-frequent perturbations from any node a fluctuation is originating from.

Power grid dynamics with and without active voltage dynamics
Some of the most fundamental dynamical network models of power grids describe a grid as a set V of N nodes (or vertices) representing synchronous machines and a set E of edges (links) representing the transmission lines between the nodes [11][12][13][14]. Jointly, the vertex and edge sets thus define a graph G = (V, E) that represents the physical structure of the grid as given by physical machines and physical transmission lines. A range of synchronous machine dynamical systems' models exist in the literature [12], with various degree of detail. They range from the second-order model that takes into account two dynamical variables per node to the sixth-order model with six dynamical variables per node. The second-order model captures the phase angle Θ n (t) of node n ∈ {1, . . . , N} and its rate of change, the node's local frequency, or phase angle velocity ω n (t) =Θ n (t) as a function of time t, both relative to a reference frame rotating at grid frequency Ω. In most Western countries, the base frequency is Ω = 2π50 Hz while some areas, including, for instance, a part of Japan have a base frequency of Ω = 2π60 Hz.
As the simplest dynamical systems' model in its class, the second-order model received broad attention in research, see, e.g. [14]. Recent advances focused on relating the physical network topology and other structural model parameters to dynamical systems responses including fluctuations, perturbations and failures [15][16][17][18][19]. A central theme is on short-term fluctuations that prevail across grids [10,[20][21][22]. Below we compare the collective dynamics of grid responses in the second-order model with temporally fixed voltage amplitudes to those in the third-order model, the lowest-order model exhibiting voltage dynamical variables. Besides our motivation to better understand similarities and differences in the collective response properties of grids with and without active voltage dynamics, we are generally interested in how the structure of the differential equation models influence non-equilibrium response properties of networked dynamical systems.
The driven second-order model of power grid dynamics is given by for all n ∈ {1, . . . , N}, where the Heaviside function is H(t) = 1 for t ⩾ 0 and H(t) = 0 otherwise, and ı is the imaginary unit. Moreover, Θ n (t) ∈ R denotes the phase angle at time t and it is first time derivative ω n (t) =Θ(t) n the local phase angle frequency. The quantity δ nk is the Kronecker delta function, indicating that only node k is externally driven. As we focus on linear response theory below, first order approximations of network responses to distributed driving signals acting on several nodes and containing various frequencies then follow by the superposition of individual responses. The remaining machine parameters are the mechanical input or output power P n ∈ R (negative for consumers and positive for generators), the inertia β n > 0, the mechanical damping α n > 0. Sinusoidal driving is acting on one machine k with amplitude ε and frequency ω (p) . The coupling strength between synchronous machine n and m is the maximal transmissible power between them. The voltage amplitudes E * n and E * m are constants fixed in time and both multiply a susceptance B nm > 0 for each transmission line connecting nodes n and m, i.e. present in the network. Note, that the model class in this study assumes lossless transmission, i.e conductance G nm = 0 between all pairs of machines n and m. The matrix B − diag(B) constitutes a weighted adjacency matrix of the graph G that represents the power grid with its nodes and transmission lines, see figure 1(a). Recent research [22] revealed three distinct response regimes that are delineated by the collection of eigenfrequencies at the normal operating state of the second-order grid model. We reiterate the main findings, see figures 1(a) and (b). At low driving frequencies below the eigenspectrum, the grid responds homogeneously with all nodes responding to fluctuations with almost the same amplitudes. In this bulk regime the amplitudes only slightly change with (sufficiently small) driving frequency. For intermediate driving frequencies within the band of eigenfrequencies of the system, pronounced resonances emerge depending also on the exact driving and response node location within the network topology and the operating state the system is perturbed from, see [10,22] for details. At high driving frequencies ω (p) beyond the system's eigenspectrum, the collective response dynamics exhibits localization [18,22]. As a consequence, nodes away from a driven node respond noticeably weaker and the response amplitudes also decay with increasing driving frequency ω (p) , see figures 1(b) and (c).
Dynamical voltage amplitudes have been proven to be crucial for system stability [23,24]. In particular the lack of sufficiently large voltage amplitudes has been identified as the root cause in real-world power system blackouts, for instance in the northeastern United States (2003) and Athens/Greece (2004) [12,25]. We are thus motivated to generalize the previous studies to the third-order model, the lowest-order model exhibiting dynamical voltage variables.
The third-order model considers dynamical voltage amplitudes E n (t) instead of temporally fixed ones E * n , thereby extending the differential equation (2.1) of the second-order model to the system of coupled ordinary differential equations. It exhibits additional machine parameters: the voltage setpoint E f and the reactance X. For future analysis, we rewrite the coupled system of first-order and second-order differential equations as a system ofÑ = 3N first-order differential equationṡ The second-order model is recovered in the limit X → 0. In this limit, the voltage equations obtain globally attractive voltage amplitude fixed points E * n = E f , effectively reducing the system (2.5) of 3N dynamical equations to a system of 2N dynamical equations after a transient.
To ensure reliable supply of electrical power, grids are operated close to a stable fixed point, here (Θ * , 0, E * ) [23,26], with Θ * ∈ R N and E * ∈ R N + , where R + denotes the set of positive real numbers, defined via (2.5) withΘ n =ω n =Ė n = 0 for all nodes n. Numerically evaluating the linear response amplitudes for the frequency response of the third-order model (figures 1(b) and (c)) reveals three response regimes as for the second-order model: a bulk regime for low, a resonance regime for intermediate and a localized regime for large perturbation frequencies ω (p) . In the localized regime, the responses decay systematically slower with perturbation frequency ω (p) than for the second-order model, as illustrated for selected synchronous machines 4-7, see figures 1(d) and (e). This finding also implies that responses to fluctuations at any given node reach farther into the grid and away from the driven node if voltage amplitudes are dynamic.
Can we explain these observations? Specifically, what is the mechanism underlying the boosted responses and associated weakened localization? More quantitatively, how do the responses scale with driving frequency and distance between driven and response nodes?

The physical grid and shortcuts via voltage variables
To obtain (approximate) analytical evidence, we first generalize the linear response theory and its evaluation as a function of graph distance and driving frequency from the second-to the third-order model. The third-order model includes products of active voltage variables in the interaction terms (2.5). As a consequence, the linearization creates three kinds of terms, two equivalent terms with one voltage variable fixed at and the other deviating from the value at the system's operating point and one additional term with both voltages fixed and the phase variables deviating from their values at the operating point.
To make progress, we generalize our framework to analyze both second-and third-order models in the same notation. So first let us denote the set ofÑ dynamic variables in a model byṼ. To simplify notation, we equivalently refer to the same set of variables simply by integer indices 1, . . . ,Ñ. For second-order dynamics N = 2N whereasÑ = 3N for third-order dynamics. In summarized form, the dynamical evolution equation (2.5) readsẋ with x(t) ∈ RÑ, the autonomous part of the right hand side (excluding external driving) is mediated by some f : RÑ → RÑ and the driving signal by a term εg(t). TheÑ = 3N components x ν (t) of x(t) represent original variables as for ν ∈ {1, . . . ,Ñ}. The components g ν (t) are given by Edges that are present in the third-order model, but not in the second-order model are marked in orange. The graph distance i.e. if the partial derivative is not identically zero, see [27]. The setsẼ andṼ form an abstract graphG(Ṽ,Ẽ) that represents the network of interacting variables. We note thatG is a directed graph because variable interactions are not symmetric, despite the symmetric coupling in physical power grids (transmission lines work in both directions) and thus an underlying undirected physical graph G. Figure 2 schematically illustrates the relation of the graphG of dynamical variables x ν to the underlying physical graph G of nodes n for a sample power grid. For the second-order model, the networkG of variables x ν in general consists of the phase angles Θ n and phase velocities ω n constituting the verticesṼ ofG (dark blue nodes in figure 2(b)); and three types of links that constitute the edgesẼ ofG. Two types of links connect each local phase variable Θ n with the local phase velocity ω n and vice versa. These links reflect the fact that in the system of differential equations, the rate of change of ω n is given by a function that contains Θ n and vice versa.
The third type of link connects from a given phase Θ m to the phase velocities ω n at all nodes neighboring in the physical graph G, i.e. if (n, m) ∈ E. Such a link exists if Θ m appears in a coupling term sin(Θ m − Θ n ) on the right hand side of the differential equation that defines the rate of change of ω n , see equation (2.1).
For our quantitative analysis of how response amplitudes decay (see next section), we consider shortest paths p n,m between two nodes m and n in the physical graph G in relation to the shortest directed pathsp ν,µ from one variable x µ to another x ν inG. In the example illustrated in figure 2(a), a shortest path p 5,1 in the physical graph G between node m = 1 (represented by a photovoltaic array) and node n = 5 (represented by a wind power station) traverses four transmission lines in the grid and is thus of length d G (5, 1) = |p 5,1 | = 4, see figure 2(b). In the second-order model, the effective distance along the shortest pathp ν,µ on the variable graphG from the phase velocity variable x µ = ω 1 at node m = 1 to x ν = ω 5 at node n = 5 is d (2nd) G (ω 5 , ω 1 ) = |p ν,µ | = 8, see also figure 2(b). As a consequence, a fluctuation signal acting upon variable ω 1 at node 1 influences at least 8 − 1 = 7 intermediate variables before affecting ω 5 , although the physical graph distance dG(ω 5 , ω 1 ) = 4 is much smaller.
In general we find for the second-order model (2.1) that the shortest path distance dG(ω ν , ω µ ) from one phase angle velocity ω m to another ω n is given in terms of the inter-node distance d G (n, m) as m). (3.5) The active voltage variables, as they for instance appear in the third-order model, shorten these inter-variable distances. In the example of figure 2(b), the additional voltage amplitude variables reduce the shortest path distance from d (2nd) In general, such shortcuts enabled by active voltage dynamics in the third-order model result in an inter-variable distance of Interestingly, further variables added to the third-order model, as they appear, e.g. in higher-order dynamical power grid models, do not further shorten the shortest path.

Quantifying the fluctuation response boost
We now analytically quantify the different localization patterns of high-frequency responses of second-and third-order dynamics by studying them in the same generalized representation of the form (3.1).
Write the asymptotic series of the response x(t) in the form as the perturbation amplitude ε → 0, see [29]. Here the generalized state vector x * denotes the vector of all dynamical variables at the system's operating point (a linearly stable fixed point of the formal equation (3.1)). Linear response theory consists of deriving and solving differential equations for the first-order coefficients a (1) (t). For notational simplicity, we drop the superscript ' (1) ' in the following, as it indicates the order of the approximation that we keep fixed at first order. Expanding f(x) about the fixed point x * in powers of ε and substituting into (4.1), we obtain a set ofÑ coupled equations where the Jacobian J has the elements for ν, µ ∈ {1, . . . ,Ñ}. Collecting all terms first order in ε yieldṡ a(t) = Ja(t) + g(t), (4.4) which is solved by the components (see appendix) with the steady response factor and the transient response Here, V and V −1 denote the similarity transform matrices under which J takes diagonal form and λ ℓ denotes the ℓth eigenvalue of J. The index κ := k + N indicates the perturbed response variable in the network of interacting response variables, on which the response amplitudes depend. The transients (4.7) decay with time because all eigenvalues λ ℓ have negative real parts due to the (linear) stability of the fixed point x * . We are interested in the steady response amplitude |A ν |, the absolute value of the complex number A ν ∈ C, at a given variable x ν . The response regime analysis of the second-order model [22] extends to the third-order dynamics via the response amplitudes |A ν |. The resonance regime is bounded by the smallest and largest imaginary part Im(λ ℓ ) of the eigenvalues λ ℓ , as one denominator of equation (4.6) is minimized if Im(λ ℓ ) = ω (p) . The homogeneous bulk regime [22] extends to the third-order model for frequency and phase-angle dynamics, as a consequence of the invariance of the model class to collective shifts of all Θ n . Such an invariance does not extend to the voltage variables, where each voltage response variable approaches a node-specific value in the limit of ω (p) → 0 (see appendix).
Here, we focus on the localized response regime as shown in figures 1(b) and (d). To retrieve a leading order asymptotic scaling for large frequencies, as ω (p) → ∞, we determine the largest exponent of ω (p) appearing in |A ν |. The following steps (explicated in full detail in the appendix) yield the leading order approximation. First we expand each term in (4.6) to obtain a common denominator for all terms. The common denominator is a polynomial in ω (p) of degreeÑ. The second step is to determine a relation between the exponents of ω (p) and the summed exponents of all eigenvalues λ ℓ that appear in every numerator. We find that the powers of ω (p) q are appearing in terms with λÑ for which we have The latter links the shortest signal propagation path length dG, via the lowest order, non-zero matrix element (J d G ) νκ , with the largest power of (ω (p) )Ñ −1−d G with a non-zero prefactor. Because the common denominator is a polynomial of ω (p) of degreeÑ, we finally find the asymptotic leading term approximation νκ | (ω (p) ) d+1 as ω (p) → ∞. (4.11) Here, we denote the distance dG(x ν , x κ ) between perturbed response variable x κ (t) and receiving response variable x ν (t) simply by d and a complicated variable-specific factor by |Φ [d] νκ |. In particular, the asymptotic scaling relation (4.11) indicates that the response amplitudes of the phase velocity variables ω n decay exponentially with distance d = dG(ω n , ω k ) to the frequency variable ω k at the driven node k and algebraically with the driving frequency ω (p) .
For simplicity we in the following denote the amplitudes in terms of the original dynamical variables (instead of in terms of the abstract indices ν), i.e. as (4.12) Applying the graph distance mappings (3.6) for the third-order model and (3.5) for the second-order model, respectively, reveals The power law exponents at response variable ω n at synchronous machine n at distance d G from the perturbed machine k is summarized in table 1.
In particular, responses are of the same order of magnitude (same exponent) for nodes that are nearest neighbors to the driven node k (at physical graph distance d G = 1); differences in exponent emerge at d G = 2 and become more and more pronounced as the physical graph distance d G grows. More quantitatively, the Table 1. Exponents in leading term approximations of the response amplitudes in the localized response regime: deviation between second-order and third-order model in leading term exponent grows linearly in d G between synchronous machines, such that difference between second-order and third-order model becomes more expressed at larger distances dG .
2nd-order m., leading exponent dG + 1 = 2dG + 1 3 5 7 9 11 . . . 3rd-order m., leading exponent dG + 1 = dG + 2 3 4 5 6 7 . . . one-term leading-order approximation indicates a relative response strength whose exponent grows linearly with d G and (4.14) Figure 3 confirms the asymptotic scaling, showing that the response strengths and more distant machines is several order of magnitude larger, when considering dynamic voltage variables in the third-order model. We reiterate that more detailed dynamical grid models, such as the fourth-, the fifth-and the sixth-order models do not further decrease the inter-variable distance, so the boosted responses do not automatically grow further. Nevertheless, there are still options for the responses to grow (or shrink). First, the prefactors |Φ (3rd) nκ | in (4.13) do depend on the model details. Second, our approximate predictions take into account shortest paths while longer paths may jointly contribute to substantially change the responses quantitatively.

Discussion
In summary, we have illustrated that the main dynamical regimes of collective phase velocity responses to fluctuating input signals-bulk responses at low driving frequencies, heterogeneous resonant responses at intermediate, and localized responses at high frequencies-constitute generic features of grids with and without active voltage dynamics, see [10,13,14,18,22]. Further, we expanded existing theories to also cover the fluctuations of voltage dynamics, which show the same characterstic three response regimes (see figure  A1). Intriguingly, however, active voltage dynamics boosts the high-frequency responses compared to grids with fixed voltage amplitudes. As a consequence, grid responses to high-frequency inputs at a given node stay less localized and have larger amplitudes at specific other nodes in the network. Subsecond, large-frequency fluctuations have been found experimentally to be overexpressed due to the non-Gaussian nature of the frequency distribution [3]. The propagation of such fluctuations throughout the grid have been investigated theoretically with focus on the frequency increment distribution and as a function of the effective inertia in the system [5].
To explicate our theory and systematically test it against direct numerical simulations we focused on the third-order network model of power grid dynamics, the lowest-order model with voltage (amplitudes) as explicit dynamic variables, and compared it to the second-order model exhibiting constant voltage amplitudes.
Active voltage amplitude dynamics (effective reactances of substantial positive value X > 0) creates shorter effective paths between pairs of grid nodes.
The response dynamics of grid with active voltage degrees of freedom thus becomes less localized with response strengths at certain response nodes in the grid, orders of magnitude larger than in systems with fixed voltage amplitudes. We emphasize that additional variables such as appearing in the fourth, fifth or sixth order model [12,14] do not further shorten the distance between grid nodes, thus indicating the same localization exponents as reported above (4.13).
Compared to the second-order model, coupling strengths in the third-order model are often more heterogeneous, thereby increasing the importance of paths that are longer than the shortest. Analyzing the impact of such longer paths should be part of future research.
Future work on fourth, fifth and sixth order model may reveal new forms of influences by the additional variable dynamics on the localized regime, as well as on the bulk and resonance regimes. Further, the integration of renewable energy sources over long distances, such as the case for offshore wind farms, requires further DC-AC mixed grids, where the power is inserted to the AC grid with DC connectors [30]. Virtual synchronous machines where AC dynamics is artificially added to the power source's DC dynamics create new classes of nodes that should supplement the modeled synchronous generators, likely qualitatively impacting the three response regimes.
Recent work has revealed that voltage responses to fluctuating input signals (with zero average) are shifted non-linearly towards smaller voltage amplitudes in the third-order model, causing voltage collapses at a specific driving amplitude, reported as a tipping point [26]. To date, like in this work, response dynamics of power grids have been mostly studied in the linear response regime [10]. Hence, genuinely nonlinear response phenomena, as the reported nonlinearly shifted voltage amplitudes, are likely to be uncovered. It thus remains as a challenge for future theoretical research to fully reveal how active voltage variables impact the consequences of short-term, sub-second fluctuations. Of at least equal importance is to uncover where the resulting fluctuation responses may induce tipping points and thus a potential loss of reliable operation. Such research seems particularly relevant on potentially existing but unknown response phenomena.
In general, we believe that research on the collective dynamics of power grids needs to particularly address (i) the more networked topology of grids upgraded with additional transmission lines, (ii) the more and more prominent high-frequency fluctuations resulting from renewable sources and consumer dynamics, (iii) the more frequently occurring high-load situations where additional short-term fluctuations can have a larger impact on grid performance, and (iv) the generically nonlinear nonequilibrium response dynamics in strongly networked systems.

Data availability statement
No new data were created or analyzed in this study. with x ∈ RÑ and f : RÑ → RÑ, perturbation components g ν (t) = δ νκ H(t) cos(ω (p) t), (A.2) and an asymptotically stable fixed point x * with initial conditions x(0) = x * . In analogy to the one-dimensional systems, we here expand the vector-valued response dynamics into a vector-valued asymptotic power series As demonstrated in the main part, the linear approximation (first order in ε) about a fixed point x * of (A.1) takes the following formȧ with the systems' Jacobian J. In the following, we drop the superscript ' (1) 'indicating the order of approximation for simplicity of notation. Equation (A.4) constitute a system of coupled, linear differential equations with a time-dependent inhomogeneity. If the Jacobian J is diagonalizable and one node is driven by one a frequency signal, this system is reducible toÑ instances of the one-dimensional first order response differential equation and Each component B ν decays as t → ∞ because all eigenvalues λ ℓ have a negative real part.  = 1 + N). The collective response dynamics exhibits a bulk, a resonant and a localized regime, depending on the driving frequency ω (p) . The bulk regime is different for the voltage amplitude variables: they still only slightly vary with (sufficiently small) changes in ω (p) yet there magnitude varies among the nodes. Vertical gray lines mark the support of eigenspectrum of the system.

Response theory in one dimensional systems
Here, we present a known derivation of the linear response function a(t), defined via the differential equationȧ with the initial condition a 1 (0) = 0. The homogeneous solution a (h) 1 (t) of the separable, linear differential equation is given by (A.14) A particular solution can be obtained by either guessing a solution of the form C p e ıω (p) t or by the method of variations of the constant, which we explicate in the following. We take the homogeneous solution and treat the constant C as being time dependent itself C → C p (t) and inserting this function on both sides of the differential equationĊ The particular solution a (p) (t) follows as The comprehensive set of solutions of the linear differential equation is then given by the superposition of homogenous and particular solution where the constant of integration C can be determined from the initial condition a(0) = 0, yielding

Asymptotic behaviour in the localized response regime
In this section we analyze the asymptotic behavior of the response factors A ν in the limit of large perturbation frequencies ω (p) ≫ max(Im(λ ℓ )). We aim to find the term with the largest power of the perturbation frequency ω (p) , which is the leading term in the localized response regime. First, we bring all eigenmodes in equation (A.11) to a common denominator ∆, which yields The common denominator ∆ is a polynomial of ω (p) of degreeÑ. Expanding the product and rewriting it as a polynomial of ω (p) yieldsÑ (A.20) The largest powers of λ ℓ occurring in C [ℓ] l relate to powers of the Jacobian and therefore to all walks between responding and perturbed variable of the same length as the power. To retrieve a leading term approximation we therefore need to express the largest power of λ ℓ in C (A.23) The first term is a sum over all eigenvalues and therefore identical for all possible C [ℓ] l . The prefactor is indeed a polynomial of λ ℓ with degree l = 1. Under the assumption the assertion hold for l we will show that the assertion follows for l + 1. where the last term is a polynomial of λ ℓ of degree l + 1 following from the assumption that the proposition holds for l. The first term is again independent of ℓ. With these findings the response factor takes up the form: We conclude that l is the largest power of λ ℓ occurring in a term (ω (p) ) 3N−1−l . Each summand l m=0 β [l] m (λ ℓ ) m (ω (p) )Ñ −1−l (A.26) is a polynomial of λ ℓ of degree l. The largest power of the perturbation frequency ω (p) corresponds to the lowest power of λ ℓ that under similarity transformation yields a nonzero contribution where the sum is the (ν, κ)-component of J l of degree l, identifies the scaling of A ν and with that of |A ν | of a response variable that is at distance d in the network of interacting variables from the perturbed variable κ (which is ω k in our setting). It scales as νκ | (ω (p) ) d+1 as ω (p) → ∞, (A.28) where |Φ [d] νκ | is determined numerically from equation (A.25).