Spatio-temporal complexity of power-grid frequency fluctuations

Power-grid systems constitute one of the most complex man-made spatially extended structures. These operate with strict operational bounds to ensure synchrony across the grid. This is particularly relevant for power-grid frequency, which operates strictly at $50\,$Hz ($60\,$Hz). Nevertheless, small fluctuations around the mean frequency are present at very short time scales $<2$ seconds and can exhibit highly complex spatio-temporal behaviour. Here we apply superstatistical data analysis techniques to measured frequency fluctuations in the Nordic Grid. We study the increment statistics and extract the relevant time scales and superstatistical distribution functions from the data. We show that different synchronous recordings of power-grid frequency have very distinct stochastic fluctuations with different types of superstatistics at different spatial locations, and with transitions from one superstatistics to another when the time lag of the increment statistics is changed.


I. INTRODUCTION
Power-grid systems represent one of the most complex and largest man-made technological structures, which are permanently active and constantly evolving. These systems supply the power needed for modern society to function. In order to exchange power between producers and consumers, the grid systems operate as coupled oscillators in strict phase locking, rotating synchronously at a nominal power-grid frequency (e.g. 50 Hz or 60 Hz). Understanding the nature of power-grid frequency is crucial as it represents one of the key observables in power system operation and stability as it measures the power balance in the grid: It increases in periods of over generation and decreases in periods of scarcity [1]. If frequency deviation exceeds a threshold dedicated control power plants are ramped up or down to restore the power balance [2]. At first glance, one could be led to believe that the powergrid frequency is identical across a power grid. This, however, is not the case, as many studies have indicated [3][4][5][6][7]. Besides oscillations in areas in a power-grid, denoted as intra-area oscillations [8], and equivalently across areas in a grid, denoted as inter-area oscillations [9], the presence of other small-scale fluctuations is ubiquitous. These stochastic fluctuations, i.e., stochastic noise with complex spatio-temporal properties, carry their own physical relevance and dictate particular physical aspects of each location's properties. They are influenced by many different factors, such as demand fluctuations, fluctuations in renewable energy production, control actions, trading, and so on.
The spatio-temporal complexity inherent in small frequency fluctuations in power grids is immense, and re-quires new techniques of analysis to obtain further insight. In this article we apply superstatistical analysis to this problem, concentrating onto the increment statistics of experimentally measured frequency fluctuations. Superstatistical methods, as introduced in Refs. [10,11], provide a general approach to describe the dynamics of complex non-equilibrium systems with well-separated time scales. These types of models generate heavytailed non-Gaussian distributions by a simple mechanism, namely the superposition of simpler distributions whose relevant parameters are random variables, fluctuating on a much larger time scale. Originating in turbulence modelling [12], superstatistics has been applied to many physical systems, such as plasma physics [13,14], spin systems [15], cosmic rays [16,17], self-gravitating systems [18], solar wind [19], high-energy scattering processes [20][21][22], ultracold gases [23] and non-Gaussian diffusion processes in biophysical systems [24][25][26]. Furthermore, the concept has been successfully applied to other areas as well, for example to environmental time series (e.g. oxygen concentration in rivers) [27], wind statistics [28], air pollution [29], bacterial DNA [30], financial time series [31,32], rainfall statistics [33], train delays [34], and plane delays [35]. In all these cases, an underlying simple distribution, typically Gaussian or exponential, is identified to explain the observed heavy tails of the marginal distributions when integrating over a fluctuating parameter. These tails often decay with a power law, although the precise form of the tails depends on the kind of superstatistics considered [36].
One of the most commonly used methods to examine the properties of time series data on smallest temporal scales is to study their increment statistics [37,38], a technique very well-known from turbulent flows [11,12,[39][40][41]. These increments-a set of differences of the time series with given fixed time lag-carry the most fundamental stochastic properties of the underlying noise within the data. For our particular application, these properties include stochastic fluctuations, their scaling properties, and correlations between fluctuations. One common aspect is the emergence of heavy tailed statistics in incremental time series, a property well suited to be described by superstatistics. With respect to powergrid frequency fluctuations, on the one hand we know that two synchronous measurements of the power-grid frequency in the same synchronous region will be basically identical at timescales > 5 ∼ 10 seconds. On the other, recent studies by us [42][43][44] indicate that at scales < 5 seconds the phase and amplitude synchronisation is not fully achieved, thus recordings of distant locations show local independent properties.
In this paper, we apply the superstatistical approach to unravel the underlying physics of synchronous power-grid frequency recordings, taken at six different sites in the Nordic Grid. We examine the characteristics of superstatistical properties in incremental time series of the powergrid frequency. We show that, as observed in Ref. [44], at scales roughly > 2 seconds the strong phase locking in coupled power-grid systems dictates the statistical properties of power-grid frequency. At small timescales, < 2 seconds, each recording, at a different location, shows distinct superstatistical properties. We also observe that below < 2 seconds there is no change of the entropic index q, suggesting no change in the superstatistical description in that range. Generally the southern region of the Nordic Grid exhibits higher entropic indices q, indicating that the internal properties of the processes vary greatly.
Our main result, from the superstatistical perspective of the complexity of the power grid, is that the grid consists of different spatial regions with different types of superstatistics, which themselves depend on the time lag chosen for the increment statistics. We observe best fits for a combination of lognormal, Gamma (or χ 2 ), inverse Gamma (or inverse χ 2 ), and F -superstatistics in various spatial regions. This complexity is higher than in, e.g. isotropic turbulent flows, where usually only one type of superstatistics is observed, which typically is lognormal superstatistics [12].

II. BACKGROUND AND METHODS OF ANALYSIS
A. Power-grid frequency dynamics In this article we focus on six synchronous power-grid frequency recordings from the Nordic Grid recorded continuously between 21:00 of the 9 th to 09:00 of the 11 th of September, 2013, with a time sampling of 0.02 seconds, sufficient to examine the stochastic properties at each location of the recordings in great detail. The locations of the recordings are indicated by the acronyms of the universities where the recordings were taken: Chalmers University of Technology Gothenburg (CTH); Faculty of Engineering, Lund University (LTH); Royal Institute of Technology Stockholm (KTH); Luleå University of Technology (LTU); Tampere University of Technology (TTY); Aalto University (AU). These data come from a former phase-measurement unit network of collaborating universities in the nordic countries, which was in operation between 2012 and 2014 [45].
In Fig. 1a we show the approximate location of each recording illustrated on a map of the Nordic Grid synchronous area (comprising Norway, Sweden, Finland, and Zealand in Denmark). In panel b we show excerpts of half an hour of recordings, vertically displaced for clarity.
Synchronous power-grid systems operate at a set frequency to ensure synchrony across the synchronous region-this is the nominal frequency of 50 Hz (60 Hz). From a physical point-of-view, we can understand each location as a node on a complex network of phase-locked oscillators, each with a given inertial mass M i , a damping constant c i , obeying (in a very reductive description) with θ i the angle of each machine in a co-rotating frame of W interacting machines, P m i the power produced and injected by that machine, and P e j the electrical power extracted by all other machines/loads. From a more practical point-of-view, we do not have access to this information for all synchronous machines i in a given area. What we can easily extract from a power-grid is, for example, the frequency f i = dθ i /dt (or angular velocity) at a particular location. This is given, in good approximation, by a Langevin equation [46] df with γ i a friction force (related in some fashion to c i and M i ) and Γ i (t) some noise function, possibly temporally and spatially correlated, with amplitude σ i .
Here we are not directly interested in the power-grid frequency, but in the incremental properties of the six measured time series. The increments of the recordings ∆f τ (t) are given by [38] ∆f where τ is the incremental time lag. In this way we move from a picture in real time t to a "scale process" in incremental time lag τ , particularly for small τ < 5 s. Note here that by examining incremental processes, on these short time scales, one basically excludes the deterministic elements of the process and deals solely with the stochastic characteristics of the fluctuations themselves. This eliminates concerns about deterministic activity or other long-scale phenomena (> 5 s), e.g. dispatch activity, control mechanisms, or power flow changes.

B. Leptokurtic increment statistics
One of the most important properties to study for such incremental time series is their probability distribution. Some previous studies of power-grid frequency involving a stochastic element already remark on the distinct features seen at the level of the increments [3-5, 43, 44]. Increments statistics often is non-Gaussian, displaying heavy tails for small time lags, commonly quantified by the kurtosis κ(∆f τ ), i.e., the fourth standardised moment, which is given by where X is a random variable (representing the increments), E[·] is the expected value, µ X is the mean value of X, and σ 2 X its variance. In Fig. 2a-c we display the distributions of the increments at the six locations at the shortest incremental lag τ = 0.02 s and at τ = 1.20 s in a vertical logarithmic scale. Firstly, we note that rather distinct distributions are immediately evident-some with considerable heavy tails, some less so. Recall that a normal distribution, in a vertical logarithmic plot, is an inverted parabola. This is a first, yet clear evidence for the presence of different local stochastic properties in power-grid frequency fluctuations. There can be numerous reasons for this: characteristics of the local generation, e.g. renewable energy generation, nuclear, or fossil fuel, differences of local consumption patterns, distant consumption requiring transferring power over transmission lines, and so on.
In order to adequately quantify the heavy-tailedness of these incremental time series we examine the kurtosis κ as a function of the incremental lag τ , in a similar way as this is done in turbulent flows [47]. In Fig. 2d we show the kurtosis κ(∆f τ ) of the incremental time series ∆f τ with τ ∈ [0.02 s, 5 s]. The horizontal dotted line indicates the kurtosis of a normal distribution, i.e., κ N = 3. One can clearly observe all recordings exhibit κ(∆f τ ) > 3 for all τ , and as τ increases, the kurtosis of the increments tends to κ(∆f τ 0 ) ≈ 3.35. Distributions with large kurtosis, i.e., κ > 3, are called leptokurtic (conversely platykurtic). Similar phenomena were seen for various other powergrid frequency recordings [42,43], yet a clear explanation why the incremental time series does not fully relax to a normal distributions with κ = 3 on a large scale is currently absent.
In the following section we will employ superstatistics as the possible mechanism to explain the presence of leptokurtic incremental time series. Such a superstatistical approach offers an explanation for the observed leptokurtic probability distributions. However, for this to work one has to carefully examine if time scale separation is realised for power-grid frequency fluctuations, such that a superstatistical treatment is justified.

C. Superstatistical generation of leptokurtic distributions
We start from the usual assumption of superstatistics that the underlying equilibrium state of the increment statistics is a simple distribution, i.e. Gaussian, on a suitable time scale T to be determined. This agrees with our formulation of the power-grid frequency dynamics as a Langevin equation in (2).
The increment statistics of the entire time series can then be viewed as a superposition of Gaussian distributions with different variances, weighted via a scaling function f (β), in itself a normalised probability distribution, such that the probability density function p(∆f τ ) of the increments ∆f τ is given by with i.e., a normal distribution dependent on a scaling parameter β, regarded as a volatility. What underlies the description of a superposition of Gaussian distribution has a clear physical explanation: The process' internal characteristics change in time. This is not hard to imagine. The amount of power generation and consumption changes over the day, and so does the contribution of each type of energy source, the total inertia of the system (linked to the number of conventional generator connected to the power grid), amongst many other properties of generation, consumption, and power transport.
Two points are crucial here: firstly, the superstatistical distribution f (β) can be different at each spatial location, as we will see. Secondly, in principle, the shape of f (β) can change for different incremental lags τ . The probability density function f (β) can generally be any appropriate and normalised distribution with support in [0, ∞). We will subsequently discuss four candidate distributions for f (β): the log-normal distribution, the Gamma distribution (or χ 2 distribution), the inverse-Gamma distribution (or inverse χ 2 distribution), and the F distribution, in line with what was done in Ref. [10]. In our analysis, we will also show that it is hard to extract a particular distribution f (β) in a unambiguous way, i.e., many distributions are consistent with the data.
The procedure to extract f (β) from a given data set is as follows. Under the assumption that different scales exist in the incremental time series of power-grid frequency recordings, one segments the data into small sub-slices, called 'snippets' in the following. Then one studies each snippet's probability distribution. Particularly, one examines each snippet of the incremental time series-for τ = 0.04 s to begin with-and extracts its kurtosis, which we then average for all snippets with a window of length with · δt the average over each snippets' length. For very short snippets the average snippet kurtosis k(δt) is smaller than 3, i.e., the average distribution shows a lack of tailedness (it is platykurtic). As the snippets increase in size δt, the average kurtosis grows. As we have seen in Fig. 2d, in the limit of the snippet being equal to the entire time series, the kurtosis is larger than 3, thus there is a certain snippet size-at the long superstatistical time T -at which the average snippet kurtosis is 3 and thus the snippets are, on average, Gaussian distributed for this particular time scale. In Fig. 3 the procedure of finding the long superstatistical time T is illustrated, for the incremental time series with τ = 0.04 s. The circles indicate the crossing at which the snippets are normally distributed, from which we determine the long superstatistical time T . Notice here that this long superstatistical time T is depending on the location where the time series is measured. It also varies for different incremental lags τ .
Having obtained the superstatistical timescale T for each time series (see Fig. 3), i.e., the snippet length at which the average snippet kurtosis is 3, one can subsequently extract the distribution of the scaling function f (β) by extracting the inverse of the variance of each snippet at δt = T We thus get a distribution of values for β T (t) from which we determine the scaling function f (β) by simply finding its distribution, i.e., examining its histogram. Notice that if all snippets have the same variance then the function f (β) is very narrow (in the limit it is a single point, i.e., a Dirac delta function). This means that there are no changes of β and the internal characteristics of the time series remains the same over time. If the variance of each snippet varies, one obtains a distribution, which is given by f (β). In Fig. 4 the underlying distributions f (β) are extracted from the incremental time series at incremental lag τ = 0.04 s [11], simply by doing a histogram of β as observed in each time slice. In panel a we display the distributions off (β) = f (β)/ max(f (β)), which have been rescaled such that each distributions' peak has a maximum value 1, so that they are visually comparable. We immediately see different widths for f (β) for the different locations.
As stated, we do not assume a priori a specific distribution f (β)-on the contrary, we wish to find this just from the distribution of the volatilities β in our given data set. Note that many theoretical distributions can be compatible with the data. We will show that many suggested distributions are compatible, although the Kolmogorov Smirnov test will point to a particular one as having least distance. One then sees that different location have different f (β) distributions.
An interesting result is the fact that the shape of the distribution f (β) is influenced by the time lag of the increment statistics. Accordingly, at larger time lags one obtains different optimal fits. This is shown in Tab. II. We consider here four candidate distributions. A lognormal distribution f logN (β) of two parameters s > 0 and µ with probability density function given by with β ∈ (0, ∞). A Gamma distribution f Γ (β) of two parameters b > 0 and c > 0 with probability density function given by with β ∈ (0, ∞) and Γ the Gamma function. An inverse Gamma distribution f invΓ (β) of two parameters b > 0 and c > 0 with probability density function given by with β ∈ [0, ∞). Lastly, we consider as well an F distribution f F (β) of three parameters, v and w positive integers  (14), and relative percentual difference from smallest Dn (lower value) for volatilities β at incremental lag τ = 1.20 s. The two locations LTH and KTH are not examined in Sec. III due to their coupling, leading to large variations in each recordings' kurtosis. Note the differences of optimal fits as compared to Tab. I, which is for a much smaller time lag.
and b > 0, with probability density function given by Naturally the subsequent question is how to assess which distribution best fits the volatilities β we extracted. The Kolmogorov-Smirnov test is a non-parametric test that evaluates the equality of two continuous distribution functions via their cumulative density functions (CDF). The Kolmogorov-Smirnov distance gives the maximal difference between the empirical cumulative distribution function F n (β) (for us given by the the CDF of the empirical f (β)) and a chosen cumulative distribution function F (β) (one of our four candidate distributions). It is given by with sup denoting the supremum. We do not have access to the true cumulative distribution function F n (β) thus we employ a numerical maximum likelihood estimation for each of the four aforementioned distributions (10), (11), (12), and (13). In Tab. I we display the Kolmogorov-Smirnov statistics D n in (14) for each location and all four distributions for the incremental lag τ = 0.04 s. We note two things: Different distributions are a best fit for different locations, but often the Kolmogorov-Smirnov distance D n is considerably small for another distribution as well. The highlighted numbers indicate the smallest Kolmogorov-Smirnov distances D n . These are compared percentagewise by taking the difference of the Kolmogorov-Smirnov distances D n divided by the smallest D n of each location. Whereas in CTH a Gamma distribution minimises the Kolmogorov-Smirnov distance D n , in LTH the Fdistribution is better suited, in KTH and LTU a Gamma distribution is again suited best and lastly in TTY and AU the inverse Gamma distribution comes first. Apart from AU, to which we will return later, there is always another distribution which would fit the empirical distribution f (β) with a comparably small Kolmogorov-Smirnov distance D n . In other words, the best theoretical model cannot be identified unambiguously.
We also note that the best-fitting distribution can change if the incremental time lag τ is changed, e.g. compare Tab. I with Tab. II. This is not surprising for two reasons: firstly, the strong coupling between the locations can influence the statistics of each location. Secondly, the data itself is limited to 36 hours of continuous measurements, which is not sufficient to uncover with exactness the underlying statistics, if there truly is a single one.
We thus conclude that that there is some ambiguity to identify the precise particular form of superstatistics from the given data. This nevertheless does not prevent us from examining in detail intrinsic properties of the incremental time series as a function of the incremental lag τ , as we will do in the next section.

III. SUPERSTATISTICAL PROPERTIES AS A FUNCTION OF THE INCREMENTAL TIME LAG
In order to ensure that a superstatistical description is viable, we need to determine whether a given time series achieves a local equilibrium at a scale much smaller than the superstatistical variation time scale of the parameter β, given by the long superstatistical time T . To do so, we examine the relaxation time of the correlation function of the incremental time series. Recall that the auto-correlation function is given by where for our case here X(t) = ∆f τ (t) (or X(t) = β(t)), i.e., we examine the auto-correlation function of the incremental time series ∆f τ (t) (or of the volatilities β(t)).
Assuming that the correlation function initially decays exponentially for the incremental time series, we can extract the decay time ρ of the exponential relaxation, i.e., ρ such that C(d) = e −1 C(0), which dictates the short superstatistical time d.
To ensure a superstatistical description is possible, the short superstatistical time needs to be smaller than the large superstatistical time d T , as the names suggest [48]. This guarantees that, locally, each incremental time series reaches an equilibrium before the larger scale superstatistics changes the physics of the process. In Fig. 5 we examine this relation for varying incremental lags τ , confirming that for incremental lags τ 1. superstatistics loses validity, with the clear exception of the recordings at AU. We note that in another work by us [44] we show that this is roughly the same scale where the power-grid frequency increments at different location lose their independence and their phases become effectively identical to each other. This seems to be simultaneously the scale where the superstatistical modelling loses validity. Interestingly, the ratio of statistical times T /d decreases with τ (with the exception of AU), which appears to be in contrast to what is observed in turbulent flows [11]. We have thus uncovered a time scale separation: at very short incremental lags τ < 1.2 s a superstatistical study is justified (d T ). At incremental lags τ > 1.2 s this clear time scale separation ceases to exist. In the following we will present our numerical results for the probability density f (β) at very short lags τ = 0.06 s, at large lags τ = 1.20 s, and for four intermediate choices indicated in Fig. 5 by circles, τ CTH = 0.30 s, τ LTU = 0.70 s, τ TTY = 0.16 s, and τ AU = 0.38 s. These are the inflection points observed in Fig. 5. In Fig. 6 we display the results. LTH and KTH are not presented since they exhibit a varying kurtosis and strong correlation, and a somewhat atypical behaviour.
As mentioned before, we cannot unambiguously distinguish between different theoretically possible f (β) for the different time lags, as for some distributions the Kolmogorov-Smirnov distance D n (or other metrics to determine the agreement of a fitting function with the data) are similar. We note that f (β) varies both with the time lag as well as with the location where the mea-surements are done. Take for example the recordings at AU: at τ = 0.04 s and at τ = 0.06 s f (β) it resembles an inverse Gamma distribution; at τ = 0.38 s f F (β) an F -distribution; at τ = 1.20 s again an inverse Gamma distribution. Assuming this is a a stable result and not just a statistical fluctuation, this means there are transitions from one superstatistics to another. This phenomenon of transitions between different types of superstatistics giving optimal fits to the data as a function of the time lag has been previously observed in Refs. [31,49] for financial time series (share price differences) as a function of the time lag.
Next, we display for each of the four examined incremental time series correlation functions, both for the original incremental time series and as well as for the extracted volatility β(t) (Fig. 6, bottom panels). From the physical point of view of superstatistics, the correlation of β(t) must be longer than that of the original incremental time series, as we have previously discussed. This guarantees that each incremental time series de-correlates faster than the changes in the superstatistical environment, ensuring an equilibrium is obtained locally before the system's physics changes. This is verified for all cases.
Let us finally quantify the strength of the fluctuations of the volatilities β. If there were no fluctuations of β, the distribution would be sharply peaked, i.e., in (5), f (β) = δ(β − β ). Following the original work [10], we quantify the width of the fluctuations of β via the general entropic index q, defined for any superstatistics as which evaluates the width of the variations of β. As mentioned before, if we observe no variations in β, our data would not exhibit heavy tails and the entropic index would be just q = 1.
In general one cannot analytically evaluate (5). Depending on the different distributions given in (10)-(13), one might or might not be able to solve the integral. Nevertheless, for small fluctuations of β the integral in (5) can be expanded. For small σ 2 = β 2 − β 2 we obtain with q as given in (16). The expectations β k in (17) are either directly formed with f (β) (type-A superstatistics) or with a slightly deformed f (β) ∼ √ βf (β) (type-B superstatistics), see Ref. [10] for more details. Again, if q = 1, p(∆f τ ) = p N (∆f τ | β ) and we then expect no heavy tails in our increment distribution.
In Tab. III the entropic indices q for CTH, LTU, TTY, and AU, for three different incremental time lags, are shown. The q-values do not change significantly with τ , which is understandable given that the correlation times  of the volatilities β are far larger than the chosen ranges of τ . What is noticeable is that each location has a very different entropic index q. This tells us that the width of fluctuations of the volatilities at each location are very different in amplitude, a spatial heterogeneity for the power grid on large spatial scales. Whereas large q are observed at CTH, indicating strong fluctuations in β, at AU these fluctuations are much smaller.

IV. CONCLUSION
In this article we focused on examining the spatiotemporal complexity of the stochastic properties of six synchronous recordings of power-grid frequency in the Nordic Grid synchronous area, for some example time series measured in the year 2013. A priori one would expect power-grid frequency to be essentially indistinguishable across a synchronous power grid due to the strong phase locking at play at each power generator in the power grid. This nevertheless does not preclude stochastic fluctuations being present in the recordings-these are indeed ubiquitous and they exhibit complex behaviour on various temporal and spatial scales.
We have shown that the increment probability densities at six spatially distant locations show heavy tails which are quite different for each of these locations. We introduced a superstatsitical model for this, considering the incremental statistics as arising from a superposition of Gaussian distributions via a superstatistical scaling function (a probability density of local variances). We showed that there is time scale separation in the system (a necessary condition for superstatistics to work) for small time lags, and that there is a phase transition at incremental lags τ = 1 ∼ 2 seconds where this description looses validity. This is in line with previous observa-tions on the emergence of phase synchronisation in the same recordings [44]. Moreover, we note that although we observe a fast decrease of the kurtosis as we increase each location's incremental time lag, the superstatistical properties remain very similar, i.e., the entropic indices q remain constant. This indicates that the decrease of the kurtosis for these recordings is not caused by a change in the internal properties of the system, it is more likely due to the strong phase locking between generator machines.
Although the data have a very high temporal resolution (0.02 seconds) they cover only 36 hours of activities. We nevertheless can quantify the strength of the fluctuations of the parameter β of the underlying physical process by studying the entropic index q of each recording, which is a measure for the typical width of the fluctuations in β. We observe that the entropic index q remains largely constant in the range < 1.2 seconds, yet it is different for each spatial location, indicating that the physics of the increments is highly inhomogeneous on a large spatial scale and dependent on the location of the recording in the power grid. This indicates that the underlying physics of the increment statistics is different at each location, and these fluctuations are stronger in the southern part of the Nordic Grid-the more populated area of the Nordic grid. This points to the need of having stronger local control mechanisms at points where more power is consumed.
From a superstatistical point of view, a particularly interesting result of our investigation is that the performed Kolmogorov-Smirnov tests point to different types of superstatistics being relevant at the six different locations, meaning there are different superstatistical distribution functions f (β) giving the best fits to the data in the different regions. Of course, this spatial inhomogeneity should still be confirmed with bigger data sets and longer time series of other power grids in the future. Nevertheless, it is apparently a phenomenon that deserves further investigation, as it is a phenomenon more complex than for homogeneous hydrodynamic turbulence, where usually just one type of superstatistics (lognormal su-perstatistics) is sufficient. Phase transitions from one superstatistics to another occur also as a function of the time lag for our data, this phenomenon has been previously observed for financial time series in Refs. [31,49]. Of course, some caution must be taken here, in the sense that a deeper statistical analysis for larger and longer data sets is necessary to confirm this result of different types of superstatistics acting in a co-existing way.
To summarise, power-grid frequency, recorded synchronously across the Nordic Grid in our investigation, exhibits very complex spatio-temporal behaviour, if the small fluctuations around the mean are carefully taken into account by doing increment statistics. On the one hand, each recording is very similar given the strong phase-locking within the power grid, on the other hand, the recorded power-grid frequency fluctuations follow different types superstatistics, depending on the spatial location, and depending on the time lag considered.
Conflicts of Interest The authors declare no conflict of interest.