{'fraction_non_alphanumeric': 0.03513194944666361, 'fraction_numerical': 0.00929867068299391, 'mean_word_length': 4.623826183023384, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Predictability and suppression of extreme events in complex systems

In many complex systems, large events are believed to follow power-law, scale-free probability distributions, so that the extreme, catastrophic events are unpredictable. Here, we study coupled chaotic oscillators that display extreme events. The mechanism responsible for the rare, largest events makes them distinct and their distribution deviates from a power-law. Based on this mechanism identification, we show that it is possible to forecast in real time an impending extreme event. Once forecasted, we also show that extreme events can be suppressed by applying tiny perturbations to the system.

Extreme events are increasingly attracting the attention of scientists and decision makers because of their impact on society [1][2][3][4], which is exacerbated by our increasing global interconnectivity. Examples of extreme events include financial crises, environmental and industrial accidents, epidemics and blackouts [5]. From a scientific view point, extreme events are interesting because they often reveal underlying, often hidden, organizing principles [6][7][8]. In turn, these organizing principles may enable forecasting and control of extreme events.
Some progress along these lines has emerged in studies of complex systems composed of many interacting entities. For example, it was found recently that complex systems with two or more stable states may undergo a bifurcation causing a transition between these states that is associated with an extreme event [9,10]. Critical slowing down and/or increased variability of measureable system quantities near the bifurcation point open up the possibility of forecasting an impending event, as observed in laboratory-replicated populations of budding yeast [11].
An open question is whether other underlying behaviors cause extreme events. One possible scenario is when the system varies in time and is organized by attracting sets in phase space. For example, a recent model of financial systems consisting of coupled, stochastically-driven, linear mappings [12] shows so-called bubbling behavior, where a bubble -an extreme event -corresponds to a large temporary excursion of the system state away from a nominal value. In this example, the event-size distribution follows a power law, having a "fat" tail that describes the significant likelihood of extreme events. One main characteristic of such distributions is that they are scalefree, which means that events of arbitrarily large sizes are caused by the same dynamical mechanisms governing the occurrence of small-and intermediate-size events, leading to an impossibility of forecasting [13][14][15][16][17].
In contrast, the new concept of "dragon-kings" (DKs) emphasizes that the most extreme events often do not belong to a scale-free distribution [13]. DKs are outliers, which possess distinct formation mechanisms [18]. Such specific underlying mechanisms open the possibility that DKs can be forecasted, allowing for suppression and control. Here, we show DK-type statistics occurring in an electronic circuit that has an underlying time-varying dynamics identified to belong to a more general class of complex systems. Moreover, we identify the mechanism leading to the DKs and show that they can be forecasted in real time, and even suppressed by the application of tiny and occasional perturbations. The mechanism responsible for DKs in this specific system is attractor bubbling. As explained below, we argue that attractor bubbling is a generic behavior appearing in networks of coupled oscillators, and that DKs and extreme events are likely in these extended systems.
The large class of spatially extended coupled oscillator networks covers the physics of earthquakes [19], biological systems such as the collective phase synchronization in brain activity [20], and even of financial systems made of interacting investors with threshold decisions and herding tendencies [21]. Many coupled-oscillator system models exhibiting chaos have invariant manifolds -subspaces of the entire phase space on which the system trajectory can reside. o another attractor, eventually followed by reinjection to the dominating attractor in many situations. In sum, attractor bubbling associated with riddled basins of attraction is a generic mechanism for DKs, which we conjecture apply to a large class of spatially extended deterministic and stochastic nonlinear systems. Such manifolds commonly occur in models where identical chaotic systems are coupled and synchronize. Furthermore, when invariant manifolds contain chaotic orbits, they can lead to attractor bubbling (as well as riddled basins and on-off intermittency) [22][23][24][25]. Attractor bubbling is a situation where the system trajectory irregularly and briefly leaves the vicinity of an invariant manifold containing a chaotic attractor as a result of an occasional noise-induced jump into a region where orbits are locally repelled from the invariant manifold. The system state then follows an orbit that moves away from the invariant manifold, but eventually returns to the attractor. These excursions of the system state to phase space regions far from the invariant manifold are our extreme events.
To highlight the connection between attractor bubbling and DKs we study two nearly identical unidirectionally-coupled chaotic electronic circuits in a master (mnemonic M ) and slave (S) configuration. The state of each circuit is described by a three-dimensional (3D) vector whose components are related to the two voltages and the current of each circuit (see Fig. S1 [26]). The temporal evolution of the state vectors is governed by the differential equationṡ where the dot over a variable means differentiation with respect to time, F [x] is the flow for each subsystem, c controls the interaction strength between the subsystems, and K is the coupling matrix. In general, the coupled system resides in a 6D phase space spanned by (x M , x S ). However, for appropriate values of c and K, the coupled oscillators synchronize their behavior [23], which corresponds to x M = x S . Hence, the coupled-system trajectory resides in a restricted 3D subspace (on an invariant manifold). In this case, it is insightful to introduce new 3D state vectors that describe the behavior on the invariant manifold x = (x M + x S ) /2 and transverse to the manifold x ⊥ = (x M − x S ) /2. Synchronization corresponds to x ⊥ = 0 and x = x M , and the basin of attraction associated with the synchronized state is riddled.
Here we study a pair of electronic circuits for which the equations (1) and (2) take the forṁ for j = M, S, respectively, δ S,j is the Kronecker delta (1 if j = S and 0 if j = M ), and The values of the parameters and other details are given in [26,27]. Equations (3)(4)(5) correspond to Eqs. (1-2) with x j = (V 1j , V 2j , I j ) T , where x T denotes the transpose of vector x, and the coupling matrix K is such that the matrix entry K m,n is 1 for m = n = 2 and 0 otherwise. As discussed above, attractor bubbling occurs when noise is present (e.g., thermal noise in the electronic com-ponents), when there is a slight parameter mismatch between the oscillators (the flows of each circuit are slightly different), or when both effects are present, which is the most likely situation in an experiment. Bubbling is indicated by long excursions of high-quality synchronization (x ⊥ close to the noise level) interspersed by brief desynchronization events where x takes on a large value -an extreme event -as shown in Fig. 1(a). We illustrate the trajectory of a typical bubbling event in Fig. 1(b), which is a projection of the 6D phase space onto a 3D space containing components of the invariant manifold and of the transverse manifold. It is seen that the trajectory remains for most of the time on the invariant manifold x , but undergoes a large excursion away fro m the invariant manifold during the bubbling event. Due to the nonlinear folding of the flow, the trajectory is reinjected to the invariant manifold after the bubble.
To reveal the existence of DKs, we collect a long time series of values of |x ⊥ |, use a peak-detecting algorithm to identify the bubbling events, and create a probability density function (PDF) for the event-sizes |x ⊥ | n , defined as the largest peak-value of |x ⊥ | within a burst. The length of the time series is large enough that the observed PDFs have reached statistical convergence and are stationary, in the sense that their shape does not change appreciably with the addition of new samples. The resulting distribution is shown in Fig. 2, where the event-sizes follow approximately a power-law distribution (straight line in log-log scale) for small to moderately large sizes (0.04 < |x ⊥ | n < 1.8) with exponent −2.0 ± 0.1. We apply a Kolmogorov-Smirnov (K-S) statistical hypothesis test to check that the distribution of event sizes follows a power law in this interval. The hypothesis of a truncated power-law is rejected for the raw data because there are small but statistically significant deviations from a straight line decorating the distribution. The hypothesis of a truncated power-law is accepted with the same value of the exponent obtained in the fit if we apply a decimation of correlated data by resampling the raw data [26]. This empirical observation is substantiated by a theoretical analysis based on the statistics of the perturbations affecting trajectories near the fixed point at the origin [26]. The analysis predicts that the exponent is −2 to leading order. Moreover, the observed desynchronization events can be rationalized as being associated with the structure of the repeller around the origin, consistent with current theory of attractor bubbling.
A substantial and significant peak in the distribution and subsequent cut-off that deviates from the power law is observed for the extremely large events (|x ⊥ | n > 2.4), which we associate with dragon-kings. Interestingly, the probability mass contained in the large peak associated to the DKs is approximately equal to the integral of the PDF that would result if the power law extended to infinity. This fact suggests that the DKs are events that would belong to a power-law distribution but had their size limited by some saturation mechanism that effectively determines a maximum size for the events in the system. The K-S hypothesis test verifies that this large peak in the PDF deviates significantly from a pure powerlaw, as expected from theory, using either the raw data or the decimated data [26]. Hence, the theory developed assuming linearization near the fixed point captures the essence of the bubbling (power law with exponent −2 and dragon-king peak of the PDF), and only fails to explain the tiny structures decorating the distribution. As discussed some time ago [25,28], a bubbling event is initiated by "hot spots" within the chaotic attractor that resides on the invariant manifold. The attractor is composed of a large (likely infinite) number of unstable sets, such as unstable fixed points, unstable periodic orbits, etc. [22]. Each of these sets has an associated local transverse Lyapunov exponent [28], which describes the tendency of a trajectory to be attracted to or repelled from the invariant manifold when it is in a neighborhood of the set. A system with attractor bubbling necessarily has a distribution of local Lyapunov exponents (see Fig. S2 [26]), where at least some are repelling (value greater than zero), even though the value of the weighted average is negative (attracting). The repelling sets correspond to the hot spots on the invariant manifold.
For the coupled oscillators studied here, it was found previously that one set in particular -the unstable, saddle-type fixed point at x = 0 -is exceedingly transversely unstable and is the underlying originator of the largest bubbles [27,29]. That is, there is a very high likelihood that a bubble will occur whenever x resides in a neighborhood of the origin for some time, and the largest events (the DKs) occur when the residence time is long and the approach is close. The large bubble event shown in Fig. 1(b) clearly originates near x = 0. This observation is at the heart of the theoretical approximation for the distribution of event sizes, where we approximate the dynamics of perturbations by linearization of the equations of motion (Eqs. (3) to (5), for j = M, S) near the fixed point.
The influence of the fixed point in the dynamics also allows us to predict the occurrence of a large event by real-time observation of x M , which is equal to x when the subsystems are synchronized, and finding the times when it approaches the origin. Figure 3 shows the temporal evolution of |x M | and |x ⊥ |, where it is seen that |x M | undergoes a sustained drop and remains below an empirically determined threshold value |x M | th preceding a large bubble (spike in |x ⊥ |), where the forecasting time is denoted by t p . A smaller threshold is associated with a larger event size and hence it can be adjusted to isolate the DKs. This description and ensuing results are confirmed by numerical integration of Eqs. (3)(4)(5), which shows excellent agreement with the experimental observations [26] and demonstrates that it is possible to forecast DKs in this relatively low-dimensional complex system. . The largest, extreme event, which is part of the dragonking distribution, is preceded by a long excursion of |xM | below a threshold of |xM | th = 0.50 whose value is determined empirically by minimizing the number of false predictions.
With this scheme to forecast DKs, we design a feedback method to suppress them based on occasional proportional feedback of tiny perturbations to the slave oscillator when |x M | < |x M | th [30]. In the presence of feedback, the temporal evolution of the slave oscillator (Eq. (2)) is modified to reaḋ where θ is the Heaviside step function, and c DK (K DK ) is the feedback strength (coupling matrix) used to suppress DKs. For the purpose of illustration, we assume that it is expensive or not convenient to keep this additional feedback coupling on all the time and thus it is only active for a brief interval when necessary. Figure 4(a) shows the temporal evolution of the system in the presence of occasional feedback. When |x M | > |x M | th , no feedback is applied and the small bubbling events are allowed to proceed. On the other hand, when |x M | < |x M | th , feedback perturbations are applied that are only 3% of the system size (defined as the maximum value of |x M | ∼ 4). Such small perturbation only causes a small change in x S , yet it has a dramatic change in x ⊥ : the large bubble is suppressed. Over a long time scale, feedback is only applied 1.5% of the time, consistent with the frequency and duration of extreme events. Thus, the total perturbation size averaged over the whole time, including the intervals when the perturbation is not active, corresponds to 0.05% of the system size. As a result of this occasional feedback, we observe that the largest events, i ncluding the DKs, are entirely suppressed, as shown in the probability density function for |x ⊥ | in Fig.  4(b). It is seen that the small-to intermediate-size bub-bles are unaffected; only the events that would have a large size in the absence of control are suppressed.
Our work addresses several important questions regarding complex systems. We answer affirmatively and conclusively that: 1) a particular simple, but nontrivial system displays DKs whose event-size distribution deviates significantly upward from a power law in the tail; 2) DKs can be predicted; and 3) this predictability can be used to occasionally and efficiently activate countermeasures that suppress or mitigate the effects of DKs. An important and immediate open question is whether it is possible to easily identify the unstable sets that are primarily responsible for causing DKs in the wide variety of complex systems that are already known to have attractor bubbling or in systems that may display bubbling but it is not yet appreciated that the behavior is of this type. While a specific method that is valid in all cases is unlikely to exist, the particular example studied here demonstrates that, with some understanding of the burst mechanism, large DK-type events may potentially be avoidable by devising small, well-chosen system perturbations. Key to addressing this problem is the development of new tools for analyzing models of complex systems or for time series analysis of natural systems that can identify burst mechanisms. We suggest that the use of this knowledge to devise appropriate control strategies is a worthy pursuit given the increasing appearance of extreme events and their impact on society.
HLDSC and MO acknowledge the financial support from the Brazilian agencies Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Financiadora de Estudos e Projetos (FINEP). DJG gratefully acknowledges the financial support of the U.S. Office of Naval Research, grant # N00014-07-1-0734 and thanks Joshua Bienfang for constructing the chaotic electronic circuits. EO and DJG gratefully acknowledge the financial support of the US ARO through grant W911NF-12-1-0101 and W911NF-12-1-0099, respectively. The authors thank one of the anonymous referees for suggestions incorporated in this Letter.
A single chaotic electronic oscillator used in our study is identical to that used by Gauthier and Bienfang in their study of attractor bubbling [27]. Here, we summarize the properties of the circuit for completeness, which consists of a negative resistor wired in series with a capacitor, which is coupled to an inductor-resistor-capacitor tank circuit through a nonlinear conductance, as shown schematically in Fig. S1. Its behavior is governed by the set of dimensionless equations (3)(4)(5)  the relation x T j = (V 1 j ,V 2 j , I j ). The chaotic electronic circuits are coupled by measuring the differences of either V 1 j or V 2 j in each circuit with an instrumentation amplifier, converting to a current, and injecting the resulting signal to the slave node above one of the capacitors.

II. LOCAL LYAPUNOV EXPONENTS
Attractor bubbling occurs when the largest transverse Lyapunov exponent λ ⊥ is negative (the synchronized state is stable in the average sense), but there exist positive local transverse Lyapunov exponents η ⊥ of invariant sets embedded in the chaotic attractor. The transverse Lyapunov exponents λ ⊥ and η ⊥ are determined from the variational equation where DF [s(t )] is the Jacobian of ponents that are of interest. It is the times when η ⊥ is positive that the trajectory s(t ) is near an unstable set embedded in the attractor and a bubbling event is likely. Figure S2 shows the probability density for η ⊥ using a short integration time of 0.005 dimensionless time units. We see that the distribution is very non-Gaussian at short times, and, most importantly, is sharply peaked around η ⊥ ∼ 0.55. Analysis shows that this large peak is associated with the unstable saddle embedded in the attractor at x ⊥ = 0. It is this unstable saddle that causes trajectories to transition between the double scrolls of the attractor.
The analysis so far assumes that the components of the master and slave electronic oscillators are identical and that there is no noise in the system. We took great care to match all components to better than 1%, and the noise level is small. In this case, it is useful to suppose the existence of an effective approximate invariant manifold, but these nonideal behaviors (i.e., noise and mismatch) are required to observe bubbling. When, at some time t 0 , the trajectory gets close to an ideal unstable set (such as the unstable saddle), we can regard the nonideal effects as acting like a transverse initial perturbation, which subsequently leads the orbit to experience exponential growth according to for t − t 0 smaller than the characteristic time scale over which the trajectory remains within a neighborhood of the unstable set whose (local transverse) stability is given by η ⊥ . This is the initiation of the bubble.

III. THEORETICAL DERIVATION OF THE DISTRIBUTION
We now present a theoretical analysis for the observed power law exponent −2, and we also give an argument to explain why the distribution is truncated at a maximum event size and ends in a peak. In order to build the PDF of event sizes, let us follow the trajectories x M (t ) of the and grow exponentially until the contribution of higher-order terms become significant in the  When the orbit is not in a bursting event, noise or parameter mismatch keeps |x ⊥ (t )| non-zero. Thus, when the chaotic dynamics places the orbit near the fixed point at some time t 0 , the linear approximation to the absolute separation between the real and ideal (no noise or mismatch) trajectories along the unstable direction is where |u ⊥ (t 0 )| is determied by the small noise and/or mismatch. In order to determine B , we neglect the components of the trajectories in the stable plane, |x ⊥ (t )| ∼ = |u ⊥ (t )|, an approxima- To find an approximate value for ∆(t sat ) there are two cases to consider: i) when |u (t 0 )| > |u ⊥ (t 0 )|, the trajectories for the master and slave systems remain relatively close to each other, so they both saturate together and ii) when |u (t 0 )| < |u ⊥ (t 0 )|, the distance between the trajectories grows faster than their center and one of either the master or the slave trajectory saturates before the other one does. In case i), the trajectories are on the same side of the stable plane. Let us assume, without loss of generality, that the master trajectory saturates first, on the positive direction, hence, Using the solution for the master unstable component, Eq. (S3), we can eliminate the time to write Similar analysis for other possibilities in case i) gives Notice that the maximum value of B in case i) corresponds to u max , when |u ⊥ (t 0 )| |u (t 0 )|.
In case ii) the trajectories are always on opposite sides of the stable plane and the distance between them has a larger absolute value than their center, hence it grows faster (with the same exponent, the growth rate is proportional to the initial size) and both trajectories saturate nearly at the same time, one of them at u max and the other at −u max , so that their final separation is B = 2u max . If saturation happens close to the edge of the 3-D chaotic attractor, 2u max is the maximum size possible for a burst. Hence, the PDF for the burst sizes in case ii) is roughly 2u max , which can be interpreted as the narrow peak at the end of the distribution observed experimentally in Fig. 2. In principle, it is also possible that once one of the trajectories saturates it will quickly be reinjected, before the other one has had the time to reach saturation. In this situation the burst size is smaller than 2u max , but large r than u max for trajectories that remain on opposite sides of the stable plane.
Using Eq. (S9) for the burst size, we find an approximate form for the distribution of bursts that are of the type described by case i) above ( |u (t 0 )| > |u ⊥ (t 0 )|). We write the probability density for a burst of size B as where δ( ) is a Dirac delta function, p u (t 0 ) u (t 0 ) and p u ⊥ (t 0 ) [u ⊥ (t 0 )] are the probability densities, respectively, for the initial center and initial separation between the trajectories of the master and slave subsystems after reinjection. Equation (S11) is equivalent to Assuming that the distribution of the center of reinjection is a smooth and slowly varying function of the distance to the plane, and recalling that the initial separation between the trajectories u ⊥ (t 0 ) is small, in Eq. (S12), such that where u ⊥ (t 0 ) is the average value of the distance between reinjected trajectories, which is a constant, and C = 2u max p u (t 0 ) (0)u ⊥ (t 0 ) is also constant. This last equation (Eq. (S16)) gives the power law with expoenent −2 observed experimentaly.

IV. NUMERICAL RESULTS
All the results shown in this article are experimental observations, for which the model discussed in the previous section is in excellent agreement. Here we show a comparison of the experiment and predictions of the model, obtained by numerical integration of Eqs. (3)(4)(5). In order to induce bubbling, we modified the parameters of the slave circuit by 1%. Figure S5 shows both numerical prediction and experimental observations of the temporal evolution of |x ⊥ |. Over a short time scale, we observe that the shapes of the individual pulses of desynchronization events are very similar in both theory and experiment, as shown in Fig. S5(a), while Fig. S5(b) shows that, over a longer time scale, the statistics of these pulses in the numerical simulation also appears to be similar to the experimental observation.
In Fig. S5 colored circles show the positions and values of the largest peaks |x ⊥ | n detected within each desynchronization burst, which we define as the event-size. To simplify the algorithm to locate those points we run a search for local maxima over a smoothed time series, where a maximum is defined as a point which is higher than the 3 previous and the 3 consecutive points. Once a maximum higher than a certain limit ( |x ⊥ | > 1) is located, any maximum appearing shortly after (7 time units) is rejected, to prevent single burst from contributing with multiple event-sizes. This procedure is applied to both experimental and numerical data to generate the sequence |x ⊥ | n .
With regards to the event-size distributions, there is great similarity between experimental observations and theoretical predictions, as shown in Fig. S6. dragon-kings, we find that precise value of the cut-off in the distribution is quite sensitive to the value of |x M | th . We therefore take this as an adjustable parameter in the model and obtain close agreement between the observations and the model when we take the threshold at 0.32 for the model.

V. TEST OF STATISTICAL CONVERGENCE
The probability density function (PDF) of the experimental event sizes was constructed by normalization of a histogram with a total of nearly one million events. These events will be referred here as the "full sample" and the number of events is the full sample size. We observed that histograms constructed with smaller subsamples, with one quarter of the total number, deviate only slightly from the histogram constructed from the full sample. An approximate functional form for the PDF of event sizes was derived in Sec. III. The hypothesis that we want to test here is that the values of the normalized histograms converge to some exact PDF which is not known from first principles with the increase of the number of points in the sample. This happens if each value of the histogram at a given bin is a random variable with finite average.
To assess this convergence quantitatively we calculate the square of the residuals, defined as

VI. STATISTICAL TEST OF DISTRIBUTION
The theoretical analysis presented in Sec. III suggests that the distribution of burst sizes in our system follows a power law with exponent −2, within a certain range of burst sizes, and ends in a statistically significant deviation upward from a power law, that we identify as due to the dragon-kings. There are a number of statistical tests in the literature that aim to verify whether data observed experimentally is consistent with a given statistical model (see supple-mental Refs. [31,32]). The idea behind all those tests is to check the hypothesis (called the null hypothesis, in the language of statistics) that the empirical data was drawn from a reference distribution, prescribed by the model. We use the method known as Kolmogorov-Smirnov (K-S) test [31,32], which is based on the maximum value of the absolute difference between the cumulative reference distribution and the cumulative empirical distribution (a normalized cumulative histogram constr ucted from the observed data). Figure 2 suggests the our data is nearly a power law for event sizes in the interval 0.04 < |x ⊥ | n < 1.80, Thus, in order to perform the K-S test, we first generate the normalized histogram of event sizes restricted to this interval. Our hypothesis is that, for the data restricted to this interval, the reference distribution is a truncated power law. Figure S8 shows the empirical distribution of the original, unrestricted data (black, the same curve shown in Fig. 2); the normalized histogram constructed using only the maxima that happen to fall within the range 0.04 < |x ⊥ | n < 1.80 (blue); and the reference distribution (red) which is a truncated power law. The integral of each of these PDFs is normalized to unit probability, hence all these curves show equal area in linear scale. We construct our empirical histograms using 500 bins with uniform sizes in a linear scale, and we have no empty bins in the resulting histograms. It is usually suggested to use logarithmic bin sizes to avoid empty bins in heavy-tailed distributions, but this does not affect our case, because of the large amount of bins and data points. Using a least-squares fit to a straight line in log-log scale, we find the exponent of the truncated power law is always within the range −1.89 to −2.14, for different choices of the endpoints of the interval of |x ⊥ | n , consistent with the theoretical value −2. We fix the endpoints at the values given and use the exponent −2.0 to define the reference distribution. To check that data drawn from the reference distribution passes the test, we produced synthetic data (false maxima) following the reference distribution and construct an "empirical" distribution from a sample of these false maxima with the same size N as the true data. The PDF of false maxima is shown in green in Fig. S8, where we see that it follows the reference distribution, but with large fluctuations that appear larger at the tail of the empirical distribution.
The one-sample K-S test uses the maximum absolute value of the vertical distance D between the empirical cumulative distributions functions (ECDF) and the reference cumulative distribution function (RCDF) as a statistic to decide whether the data was drawn from the reference distribution, thus, we proceed by integrating the empirical and reference PDFs up to size |x ⊥ | n to generate the corresponding empirical CDF (ECDF) and reference CDF (RCDF), shown